A Novel Rotation-Based Standardless Calibration and Characterization Technique for Free-Space Measurements of Dielectric Material

This article presents a novel transmission-only calibration technique for free-space quasi-optical material characterization, based on rotating the sample around its axis to vary the angle of incidence under which the sample is illuminated. In contrast to common time domain approaches, each frequency point is evaluated individually. Thus, no minimum bandwidth is required and artifacts due to time gating are prevented. In this article, two methods are presented: the first is based on self-calibration, such that all error terms are obtained by the measured sample itself. The second one, which is tailored for thin samples, requires two known standards. Since plane-wave illumination cannot be assumed for highly-focused beams, an analytical model for the coupling of arbitrary paraxial beams is developed, accounting for the lateral beam shift in case of angled samples. Thus, the presented methods are not restricted to free-space beams with high Gaussicity, allowing to employ a variety of feed antennas. Measurements in the frequency range from 220 GHz to 330 GHz of a well-known alumina sample verify the different calibration methods.


I. INTRODUCTION
In the last decades, the metrological characterization of dielectric materials has been an ongoing task with high demand in the field of high-frequency engineering, with applications in automotive, aeronautics, space technology, biomedicine and many more.Accordingly, a variety of characterization methods exist to determine the complex permittivity as a fundamental physical property [1], [2], [3].Resonator-based approaches [4], for example, are particularly suitable for media with very low losses, but are limited to a few discrete frequency points and are difficult to implement in the THz region due to the small wavelength and thus required mechanical precision.Quasi-optical free-space methods, on the other hand, allow broadband, non-destructive measurements [5].In this context, they place low demands on the preparation of the sample as only a sheet of homogeneous material is needed.The required size of the sample is usually given by the beam diameter in the free-space system and can be reduced by generating a strongly focused Gaussian beam [6].At the same time, diffraction effects at the edge of the sample are reduced [7].
Due to a variety of unavoidable systematic errors in a typical free-space measurement setup, such as transmission loss and mismatch in the signal path, calibration is essential to determine accurate scattering parameters from which the permittivity of the material under test (MUT) can be extracted.
In the past, a variety of calibration methods suitable for quasi-optical systems have been presented.Examples of these are the Thru-Reflect-Line (TRL) [8], [9], Thru-Reflect-Match (TRM) [10], [11] and Gated-Reflect-Line (GRL) [12], [13] methods.Common, however, is the need for known calibration standards and exact positioning of both the MUT and the standards.For a large number of measurement scenarios, the exchange of calibration standards is impractical, such as for high-temperature measurements [14] or measurements under protective atmosphere, necessary in many THz measurements, when the effects of air humidity have to be suppressed [15].Furthermore, in the presence of system drift, the necessary recalibration during long-term measurements poses repeatability problems if the procedure requires replacing the MUT with calibration standards.Methods such as the TRL method usually require movement of the antenna positions resulting in an altered alignment of the propagating beam, which must be compensated for high precision measurements.This alignment problem is overcome in the TRM method, but it still requires a fully directional measurement setup.On the other hand, the reliability and accuracy of the time-domain gating in GRL methods, which is often used to suppress the mismatch error, depends on the uncertainties associated with filtering and windowing [16].For example, the available bandwidth must be sufficient to resolve multiple reflections between the MUT and the system due to mismatch at the lenses or antennas.In this case an evaluation of each individual frequency-dependent scattering point might become more useful.
In this work, a calibration method is presented that utilizes the signal difference between several transmission measurements at different incidence angles of the paraxial beam to obtain calibrated extraction data, requiring no other calibration material or device than a rotation stage.In contrast to most existing calibration methods, no reflection measurement is needed for the calibration procedure.Subsequently, a calibration method, which requires two known standards, is presented to allow accurate measurements of electrically thin samples.In contrast to previous work, which performs material extraction under oblique incidence [17], [18], [19], [20], the results are fully calibrated data and no time-domain gating is needed to resolve mismatch errors.
Another important aspect for high-precision material characterization, relevant to highly-focused quasi-optical systems, is the occurrence of focusing errors due to insertion of the MUT in the propagation path [21].Due to the incidence of the Gaussian-like beam at an oblique angle on the MUT, a lateral shift of the beam inevitably occurs.This problem can be described for Gaussian beams with a coupling model as known from [22], but requires horn antennas with high Gaussicity to be employed, since only the fundamental mode is considered in the model.Therefore, a more general approach is proposed here, which retains its validity for all paraxial beams and therefore allows the use of arbitrary antennas.The resulting correction terms can then be used to extract material parameters from the measured data.This article is structured as follows.Section II describes the system setup of a typical quasi-optical instrument and its associated error model.In Section III, prerequisites, that must be met to perform a good calibration, are listed.The S-parameters of a wave obliquely incident on a dielectric slab are derived in Section IV, followed by an extension of the plane wave model for Gaussian beams in Section V.The self-calibration procedure without any known standard is then presented in Section VI, whereas the calibration procedure incorporating at least two different standards is described in Section VII.Section VIII presents fully-calibrated characterization data of a well-known material sample and compares the different calibration approaches introduced in this article.Section IX concludes this article.

II. QUASI-OPTICAL SETUP AND ERROR MODEL
The free-space setup, as shown in Fig. 1, consists of two horn antennas connected to a vector network analyzer (VNA) for transmitting (Tx) and receiving (Rx) the measuring signal and a total of 4 polytetrafluoroethylene (PTFE) lenses.To generate a focused beam, lenses directly in front of the horn antennas collimate the radiated wave to a parallel beam, while the inner lenses focus the parallel beam to a spot.Note that it is not necessary for the system to produce a perfect fundamental Gaussian beam mode at the focal plane, as long as the paraxial approximation of geometrical optics is valid.Thus high-Gaussicity horn antennas are not strictly required as a source.The MUT, which is required to be a plate with two parallel faces, is mounted rotatably in the focal plane to allow angular adjustment.The reduced beam diameter at the focal plane prevents edge diffraction at the sample and adds the benefit of allowing fully-calibrated measurements of small samples in diameter.Since the MUT is introduced into the beam path at an angle, placement of absorbers or providing enough free space in the area, which is illuminated by the reflected signal, is essential.This prevents back-radiation of waves reflected by the MUT into the system and is an important prerequisite for successful calibration.Although the system described in Fig. 1 has been proven to be practical in many applications, the following calibration methods are not strictly limited to it; in fact a wide variety of focusing transmission-based systems are viable.
Typical free-space measurement setups can be described by a 12-term error model, as is the case for any two-port VNA measurement [23], [24].By placing an angled MUT inside the signal path and terminating the reflection at the MUT interface with an appropriate match, this error model is reduced to six terms, as depicted in Fig. 2; or four terms remain respectively if only the forward path is considered.These error terms must be determined within the calibration procedure.The actual transmitted signal X F can be unknown, but must be invariant over all measurements, thus X F = const..It must be ensured that possible drift of the transmitted signal, for example due to time and temperature, can be neglected.The received signal Y F in forward direction follows from the signal flow graph in  Fig. 2 according to Mason's rule [25]: with the error terms for transmission tracking E TF , the source match at the transmitter E SF and the load match at the receiver E LF .The isolation error E IF describes the direct port coupling.Note that the source match of the forward path is equal to the load match of the backwards path and vice versa, i.e.E SF = E LB and E SB = E LF applies.This error model retains its validity when the MUT is illuminated at different angles of incidence θ , with θ = 0 • being equal to normal incidence.
In case of a single active port, an invariant source and assumed reciprocity of the MUT, the following simplifications are valid: The tracking error E XTF and the isolation error E XIF differ from their counterparts E TF and E IF by incorporating the invariant source X F .In the following calibration procedure, the load match E LF and source match E SF appear only as product, so they are summarized to the source-load match E SL .If the isolation term E IF is sufficiently small, the error model reduces to a two-term problem:

III. PREREQUISITES
The reduced error term model in (2) is valid if the reflection of the sample is well matched and no transmitted signal is bypassing the sample, i.e. the isolation is high.For verification of this condition, a method for estimating the isolation and reflection match error terms based on the uncalibrated system is developed in the following.These estimates are calculated for various angles of incidence to find the angular range in which the error terms do not exceed a predefined upper limit.The procedures described in this section evaluate the mechanical setup itself; they are independent of the calibration procedure and do not have to be repeated for each measurement cycle.

A. ISOLATION
Omitting the isolation error terms E IF and E IB is valid if their contribution to the received signal is negligibly small.In case of a low-loss, low-reflection system, which is generally a sufficient approximation for free-space measurement systems, the isolation error of the system is estimated in the following.
A metallic plate of same dimensions as the MUT is placed inside the optical path at the angle θ , such that S 21 (θ ) becomes zero, i.e. transmission is blocked, obtaining the isolation measurement Y 0F .Consequently (1) reduces to Next, a thru measurement Y 1F is performed (S 21 = 1) by removing the metal plate: E TF and E SL are unknown at this point.Due to path loss, the absolute value of the transmission error E TF is always smaller than one, allowing to assume a worst-case estimate of As a result of this assumption, only the phase component of the transmission error remains, so that E TF ≈ e jϕ with unknown phase ϕ applies.By inserting (3) it follows: Within a low-reflection system, E SL can be omitted (E SL ≈ 0).Rearranging and taking the absolute value yields This condition must be satisfied for all angles of incidence θ , since the isolation changes with different orientation of the MUT with respect to the propagation direction.In (6) only the forward path is considered, but since the setup is symmetric, the result is equally valid for the reverse path.
The measured result of the isolation estimation of ( 6) for the presented setup is given in Fig. 3 as a maximum over the whole frequency range, i.e.
A minimum isolation E IF,est (θ ) of greater than 48 dB over the whole frequency band for oblique incidences of θ = 30 • to 60 • is observed.Noteworthy is the increase in amplitude for both low and high angles of incidence.For high incidence angles, this can be explained by a lower effective MUT diameter with respect to the beam diameter, resulting in increased spillover.At low angles of incidence, the power is essentially reflected back into the optical path, with some of the power scattering off the edges of the lenses and thus past the MUT.

B. REFLECTION MATCH
Another important prerequisite for the presented calibration procedure is the complete absorption of the signal after it is reflected at the MUT.For measurements at high frequencies this can be done by providing free space perpendicular to the optical path, i.e. in the reflection path, so that the free-space attenuation becomes large.For low frequencies or compact measurement setups however, reflections back into the optical path can become quite significant.In the following, a method for estimating these reflections in terms of the match error E M is described.Note that this method requires 1-port measurements and is not suited for transmission-only systems, where the reflection match must be estimated by other means.By placing a metallic plate in place of the MUT at the angle θ , the error model of the resulting 1-port signal path can be desribed as and is depicted in Fig. 4. For a metallic plate under normal incidence, with θ = 0 • , E M = 1 and under the condition of a low-reflection system, the following assumptions can be made: The product E M E S of the match error and the source match error E S can be neglected, so that 1 1−E M E S ≈ 1 and the directivity error E D is much smaller than the reflection tracking E T .This results in: By placing the metallic plate at two different, nonperpendicular angles of incidence, i.e. θ = 0 • and the offset angle θ 0 = 0 • , the match error E M is different for these two angles, but can be described in dependence of each other with E M (θ ) = E M (θ 0 ) + ΔE M .Thus the received signal can be written as By subtracting the two measurements and rewriting, the match difference ΔE M can be calculated: In the presented setup, the signal reflected from the angled metallic plate propagates some distance d abs before impinging on a arbitrary structure, typically an absorber plate.By changing the angle of incidence θ of the metallic reflector, this distance changes by If ΔE M (θ, θ 0 ) is now calculated with respect to multiple θ 0 , which are within a narrow angular range Δθ around θ , thus θ 0 − θ ≤ Δθ , an estimation of E M (θ ) can be given.The phase difference between the two reflections from θ and θ 0 is typically greater than 2π over the spectral range c 0 /d diff , as well as highly randomized due to the multiple, changing d diff .This leads to constructive and destructive interference of Y (θ ) and Y (θ 0 ), depending on the specific frequency point and d diff .
Consequently, an upper limit estimate E M,est (θ ) for the match error E M (θ ) is found with by taking the maximum of |ΔE M (θ, θ 0 )| over the whole frequency range and angular range Δθ .This procedure also minimizes estimation errors of E M,est (θ ) due to the highly spatial randomized reflection coefficient of a distributed structure like e.g.pyramid absorbers, as the worst possible case is taken as the limit E M,est .
For the presented setup, the measured upper limit is shown in Fig. 3. To obtain a good estimate, the angular range is set to Δθ = 5 • with a distance d abs ≈ 70 cm to the nearest obstacle, in this case absorber plates.For an upper limit of −40 dB the usable angular range is observed to be between θ =30 • and 60 • .At higher oblique angles of incidence, reflections back in the system increase mainly due to edge diffraction, as due to the small sample size its edges get increasingly illuminated.At lower angles of incidence, reflections from the frame and posts holding the lenses contribute significantly.Additionally, parts of the reflected beam are collected by the focusing lenses, i.e. lens 2 and 3 of Fig. 1, due to the spatial distribution of the beam.For θ = 0 • this contribution is maximized, with all of the reflected power being radiated back to the transmitter.

C. DETERMINING THE ZERO ANGULAR POSITION
Another important aspect is determining the angle of incidence of the beam on the MUT.While a small absolute angular error is negligible at low angles of incidence, i.e. θ ≈ 0 • , the relative error in transmission increases strongly at highly oblique angles of incidence.A manual alignment is therefore not sufficient and, to overcome this problem, an automated procedure is implemented.
By inserting a dielectric plate with known thickness and permittivity in the transmission path, the frequency response can be used to determine the zero angular position.However, if only a single frequency point is used to determine the transmission maximum, the results could be misleading due to the unknown error terms of the system, as these cannot be determined prior to the determination of the zero angular position.A more robust method would be to take the sum over the entire measured frequency band, but due to the  multi-reflection within the material, frequency regions exist in which the transmission increases with decreasing angles of incidence, as well as the opposite is the case.Since the material is known, the corresponding frequency regions, shown in Fig. 5, can be determined via the analytical calculation from (24), so that a function A(θ ) can be defined, which exhibits a maximum under the zero angular position: θ step is the discrete angle difference between two subsequent measurements and θ max is the largest angle of incidence considered.This method requires a rough alignment of the dielectric, especially to avoid an excessive high θ max , which can be done manually.Special attention must be paid to the measurement parameters, since the contrast in transmission at small angles is vanishing and thus the measurement setup is sensitive to noise.A measurement of A(θ ) is shown in Fig. 6 together with a polynomial fit.The zero angular position corresponds to the maximum of the fitted function.

IV. S-PARAMETERS OF A WAVE OBLIQUELY INCIDENT ON A DIELECTRIC SLAB
To predict the S-parameters of a wave obliquely incident on a known dielectric slab, a plane wave model according to [26] is employed.This provides the basis for a model of the Gaussian beam propagation (see Section V) present in the focal point of the measurement system and shall be briefly discussed at this point.Assuming plane-parallel interfaces at the dielectric slab, the reflection R and transmission T between two interfaces can be described using the Fresnel formulas [27].Without loss of generality, TE-polarization of the wave obliquely incident on the slab at the angle θ is assumed in the following.With relative permittivity ε r,M and its correlation to the refraction index n M = √ ε r,M of a non-magnetic, homogeneous, isotropic MUT and vacuum as background medium, they can be written as As shown in Fig. 7, the two interfaces lead to multiple reflections inside the MUT of thickness d.By summing up all individual Fabry-Pérot echoes, the transmitted signal S 21 (θ ) is given by with the wavenumber of free space k 0 and the wavenumber k M = ω √ ε 0 ε r,M of the wave propagting inside the MUT.The line lengths are given by with the refraction angle θ t calculated by Snell's law sin θ t = sin θ n M [27].The reference plane corresponds in this case to the entry and exit points of the wave at the MUT interfaces.In order to enable comparability between the S-parameters at different angles of incidence, a de-embedding of the reference plane to a plane independent of the angle of incidence must be applied.The two separated reference planes are thus contracted to one plane, resulting in

V. EXTENSION OF THE PLANE-WAVE MODEL FOR HIGHLY-FOCUSED GAUSSIAN BEAMS
In case of small beam diameters in relation to the thickness of the MUT, the plane-wave model described in Section IV is no longer valid [21].A more elaborate model, expanding the plane-wave model (24), is required to describe propagation inside of the MUT, and to account for the lateral beam displacement occurring when a sample is illuminated by an obliquely incident beam.This is done by weighting the coefficients t k from (24) with terms representing the amplitude distribution of the field in the focal zone.This has already been demonstrated for a wave configuration resembling the fundamental Gaussian mode [22].For arbitrary field distributions of collimated beams, e.g.generated by antennas with low Gaussicity, the propagating field in the focal spot can be decomposed into a set of orthogonal Hermite-Gaussian beam modes.The electrical field of the m, nth-Hermite-Gaussian mode E m,n (x, x 0 ) in rectangular coordinates is given by [28] E m,n (x, at the distance x = (x, y, z) from the beam waist origin x 0 = (x 0 , y 0 , z 0 ) , propagating in z-direction.In contrast to the notation in [28], the linear phase shift term exp(−jk 0 (z − z 0 )) is omitted, as E m,n (x, x 0 ) only describes the expansion of the plane-wave model for Hermite-Gaussian beams and not the wave propagation itself.H m and H n denote the Hermite polynomial of order m and n, respectively [28].The beam radius w x,y (z), the radius of the curvature R x,y (z), the Gouy phase ψ x,y (z) and the Rayleigh range z R read Besides the field strengths E m,n,0 of each mode order (m, n), this leaves five parameters to fully characterize each Hermite-Gaussian mode: the beam waist radius in x-and y-direction, w x,0 and w y,0 , respectively, and the three coordinates (x 0 , y 0 , z 0 ) of the focal point.Since these parameters are system-dependent and thus invariant for each measurement run, they must be known prior to the calibration procedure.The total propagating field is then given as the sum of all orthogonal modes [29], i.e.
To put the propagating field in context with S-parameters, the field strengths E m,n,0 have to be normalized, so that holds.
If a Gaussian-like beam is obliquely incident on a dielectric plate, multireflections occur due to the boundary conditions at the interface, see Fig. 8.The displaced beam position in x-direction of the kth transmitted beam is described by x k , while z k corresponds to the beam waist shift in propagation direction [30] in the medium M: Since the dielectric plate is rotated around the y-axis, there is no displacement in y-direction.The transmitted field E Tx is the sum of all k echoes originating from x k , where the coefficients t k of the individual beams follow from the plane-wave model (24), thus For the following considerations and without loss of generality, the field E Rx (x, x 0 ) associated with the receiver optics is aligned to the system origin and propagating in free space.It follows: The receive field E Rx and transmit field E Tx , corresponding to the right side (Rx) and left side (Tx) in Fig. 1, thus originate from a beam with identical Gaussian beam parameters, i.e.E sum .This is valid if the system is sufficiently symmetric.Otherwise, the two fields must be measured separately and decomposed into their Hermite-Gaussian modes, with independent mode parameters w x,y,0 , E m,n,0 and beam origins x 0 .
The resulting transmit S-parameter S 21 is equivalent to the field coupling, described by the integral of the Tx-field and the Rx-field in any arbitrary plane with z = const.[28], so Under perfect alignment and with no dielectric plate inserted, this integral is equal to one.Generally, this model correctly describes the amplitude and phase variations due to the beam shift and beam broadening introduced by an angled dielectric, compared to the plane-wave model.A closed-form expression of the integral in (37) can be derived according to [28], if the fundamental Gaussian mode TEM 00 is exclusively present in the system, greatly reducing the computational effort.

VI. CALIBRATION PROCEDURE WITHOUT ANY KNOWN STANDARDS
The model presented in Section IV, with the extension for highly-focused beams in Section V, describes accurately the input-output relationship of a wave impinging on a known dielectric plate in terms of an angle-dependent transmission coefficient.This describes the so-called forward model.For material characterization, however, it is of interest to find the inverse of this problem, so that an unknown complex permittivity can be extracted from measurement data.This inversion can only be solved analytically in certain cases.For example if both the reflection and transmission coefficient under normal incidence are known, the well known Nicolson-Ross-Weir (NRW) method [31], [32] or for low-loss materials [33] can be employed.Since with the presented method only transmission data are available, iterative methods are applied instead of providing an analytical solution.

A. UNIVARIATE MINIMIZATION METHOD
To obtain permittivity data with a calibrated measurement setup, (2) has to be solved for the angle-of-incidence-invariant error terms E XTF , E SL and the angle-dependent S 21 (θ ).Subsequently the complex permittivity ε r can be calculated from the angle dependent S-parameters by inversion of ( 37) with the a-priori-known physical dimensions of the MUT under the assumption of a homogeneous, isotropic material and the characterized instrument specifications in terms of the Gaussian beam parameters.This leaves three unknown parameters to be determined, requiring a minimum of three linear independent equations.Hence, three measurements at three different angles of incidence θ 1,2,3 are performed.Solving the resulting equation system yields the cost function f err (ε r ), where M(θ ) denotes the measured receive signal, which is modeled by Y (θ ) in (2): By minimizing f err (ε r ) using a suitable least-square algorithm, such as the trust-region-reflective algorithm, the complex permittivity ε r can be determined.Due to the iterative process, a good start value is needed, which can be obtained by an educated guess.In this case, measurement data under normal incidence, normalized with a thru measurement, can be sufficient for the initial guess.
The specific choice of the three angles under which the measurements are recorded is critical.In the presence of noisy measurement data, it is especially important to have a high contrast in the measurement data in order to obtain good results.Therefore, a high angular difference of the three measurements is advantageous, exploiting the usable range determined in Section III.

B. MULTIVARIATE MINIMIZATION METHOD
The univariate minimization method, as presented in the previous section, allows the measurement problem to be reduced to minimizing a function with a single variable.However, the accuracy can be improved by including additional measurements taken at further angular positions, as it reduces the influence of randomly-distributed measurement noise, albeit with the drawback of increased measurement time.Another advantage is the possibility to use not only the forward but also the reverse path for calibration and measurement, provided that a complete two-port system is available.By measuring both paths, M F and M B , it is possible to mitigate phase drift over time, as this drift usually occurs differentially in the measured signals.
A solution for the complex permittivity data is then given by minimizing all k cost functions for the forward path f err,F,k and backward path f err,B,k with Since the source signal is included in the tracking error E XTF , E XTB , they are different for both paths, while the source-load match E SL is identical in both cases, as can be seen in Fig. 2.However, it should be noted that the computational cost is significantly higher due to the necessary multivariate optimization.

VII. IMPROVING ACCURACY FOR THIN SAMPLES BY INCORPORATING KNOWN STANDARDS
Although the algorithm presented in Section VI allows for the extraction of permittivity data, its applicability strongly depends on the sample under investigation.For example a 0.2 mm thick SiO 2 sample, exhibits a maximum contrast of less than 3.4 dB in amplitude and 22 • in phase between 60 • and 30 • in the frequency range from 220 GHz to 330 GHz.
In the presence of noise, this limits the achievable accuracy, especially when the amplitude and phase contrast drops further for very thin samples.In this context a thin sample denotes a low physical thickness in comparison to the respective wavelength and a low dielectric constant.To improve the calibration accuracy in case of thin samples, known standards can be added to the measurement procedure to better estimate the error coefficients and consequently the permittivity.
In this procedure, calibrated transmission parameters are retrieved as an intermediate step; hence the subsequent step of converting the S-parameter to permittivity data, as in case of dielectric plates, is not strictly necessary.In fact, any arbitrary DUT, can be characterized in terms of its transmission coefficient, as long as it can be guaranteed that the reflection coefficient is zero or the reflection at the DUT is radiated out of the system.

A. UNIVARIATE MINIMIZATION METHOD
As is the case for the univariate minimization method calibration without any known standard, three measurement runs are sufficient: Namely, two reference measurements for the estimation of the two error coefficients E XTF , E SL and one measurement with the MUT.First, a thru measurement M 0,uncorrected is performed by removing the sample.As the phase centers of the horn antennas shift over frequency, perfect alignment in terms of the z-position for all frequency points cannot be achieved.For correction, ( 37) is scaled with The placement of a reference material with known thickness and permittivity ε r,ref at the angle of incidence θ 1 represents the second measurement M(θ 1 ), allowing to solve (2).Therefore, for the tracking error E XTF and the source-load match error E SL the following applies: By inversion of (2) the transmission S-parameter of the MUT can be obtained from the measured signal M(θ 2 ): θ 2 can be chosen arbitrarily and can therefore be equal to θ 1 , but must lie within the boundaries found in Section III.As (43) has two solutions, the correct sign for meaningful S-parameters has to be chosen.For passive, lossy media it follows that |S 21 (θ 2 )| ≤ 1, so the corresponding sign has to satisfy this condition.The complex permittivity can then be obtained by minimizing the cost function where S 21 (θ 2 ) is obtained via (43) and S 21 (θ 2 ) is equal to (37) under θ 2 .

B. MULTIVARIATE MINIMIZATION METHOD
Similar to the extension (Section VI and VI-B) of the univariate minimization method method (Section VI and VI-A), the calibration procedure with standards can also be improved by inclusion of multiple known standards.Likewise, the forward and backward path can be incorporated into the calibration procedure to mitigate differential drift of the measurement system.First, E XTF and E SL are determined by performing a thru measurement followed by measurements of a known dielectric under multiple angles of incidence θ k .This yields the following error functions: With the now-known error terms, ε r of the MUT can be determined by minimizing the following cost functions:

VIII. EXPERIMENTAL VERIFICATION
To verify the presented calibration methods, a known material sample is characterized.By calculation of the root-meansquare error (RMSE) of the calibrated measurement data to the material data obtained by an independent system, the performance of each presented calibration method is evaluated.The prerequisites presented in Section III apply for the following measurement setup.

A. EXPERIMENTAL SETUP
The measurement setup corresponds to the free-space setup described in Section II and is shown in Fig. 1.The two rectangular horn antennas are connected to a VNA via frequency-converter modules, enabling full 2-port measurement from 220 GHz to 330 GHz.This frequency range allows the use of 50.8 mm sized PTFE lenses with 75 mm focus length due to the small aperture size of the source.To share a common optical axis, all parts are mounted on an optical rail system, whereas the frequency-converter modules are fixed to an adjustable mount.This allows for precise alignment by manual search of the maximum transmission over the whole frequency range through a pinhole in place of the MUT.In addition, a motor-driven rotary precision stage allows the MUT to rotate freely 360 • around its vertical axis.With a beam diameter of less than 3.2 mm, the illuminated spot is found to be about 15 times smaller than a sample with a diameter of 50.8 mm at the lowest considered frequency.At an angle of 60 • , this ratio is reduced to about 7 with an effective size of the MUT of 25.4 mm.This ratio is sufficient to achieve good isolation, as shown in Fig. 3. On the other hand, the beam diameter at the highest frequency, i.e. the smallest expected beam size, is 2.4 mm and therefore greater than the lower limit of w 0 /λ > 0.9 for the paraxial assumption [28].

B. GAUSSIAN-HERMITE BEAM PROFILE
Prior to the calibration and measurement procedure, the setup must be characterized with regard to its Gaussian beam profile in the focal zone.Since it cannot be assumed that the beam has high Gaussicity, unless special horn antennas are used, this requires a characterization by means of a mode decomposition to determine the beam width w x,y,0 and the normalized field strengths E m,n,0 of each associated mode.The Gaussicity of some horn antennas has already been investigated [34], but due to alignment uncertainties in the measurement system and the insertion of the lenses in the beam path, especially the higher order modes are significantly altered.Consequently, the characterization must be done separately for each frequency point, as the radiation pattern of the antennas and the beam propagation changes with frequency.
In the following, only the notation for one dimension x is given, since the total field can be separated into orthogonal directions.Yet a measurement must be done for both axes, to obtain all independent coefficients.The beam waist radius w x,0 and the mode coefficients E m,0 in x-direction can be retrieved from the electrical field E (x) at the position z = 0 according to (25), which is given by In case of a symmetric measurement setup, the beam profile can be determined using a thin metal blade, as shown in Fig. 9.This method is known as the knife-edge method [35] and allows beam characterization to be performed with the measurement system fully set up, without requiring additional measurement equipment or a system modification that could alter the beam alignment.To not obstruct the area between the sample and the absorber, the metal blade setup was dismounted for the measurements of the MUT.
By moving an obstacle perpendicular to the propagation path, the amplitude profile A(x) in one dimension x, as shown in Fig. 10(a), follows the integral of (51) and can be written as: with the measured amplitude A 0 , when the metal blade is fully retracted.The multiplication of the electrical field with its complex conjugate, is done to attain meaningful mode  coefficients with respect to the Rx-and Tx-field, as introduced in Section V.The Gaussian beam profile is obtained by differentiating the measured amplitude profile and can be seen in Fig. 10(b) for 275 GHz.Due to this differentiation, measurement noise is amplified and strongly visible in the Gaussian beam profile.Thus, a curve fit, seen in Fig. 11(a), to extract the beam width w x,0 and Gaussian mode parameters E m,0 is performed on the measured amplitude profile A(x).This has in particular the advantage that the high frequency limit A x,0 gets estimated correctly, which cannot always be guaranteed for a curve fit to the Gaussian field seen in Fig. 10(b).This fit is done separately for each frequency point, as shown in Fig. 11(a) for the beam width w x,0 in the x-axis.Note that this fit also includes the normalized field strengths E m,n,0 of each considered mode, allowing the fit to follow the indentations around the main lobe seen in Fig. 10.To prevent overfitting to noisy measurement data, the number of modes considered is limited to 4. As it can be observed, the beam diameter becomes smaller at higher frequencies.A polynomial fit of order 3 approximates the extracted beam waist diameter and the mode coefficients to accommodate for the changed radiation pattern at different frequencies.Moreover it is used as input for the Gaussian extension of the plane wave model, described in Section V. Another important aspect is the location of the beam waist z 0 .Due to the variation of the phase center of horn antennas over their respective bandwidth, the parameter z 0 must be determined for each frequency point separately and cannot be assumed to be zero.As presented in [22], z 0 can be found by mechanical shifting the position of the horn antennas and lenses along the z-axis.z 0 is then found at the position of maximum transmission, as can be seen in Fig. 11(b).A polynomial fit of order 3 is again used to minimize systematic errors introduced by the uncalibrated setup.

C. MEASUREMENT RESULTS
In order to experimentally verify the presented calibration methods, a high-quality alumina (Al 2 O 3 ) sample with known dielectric properties from a prior measurement with a terahertz time-domain spectrometer (THz-TDS) [36] is characterized.The thickness d = 0.625 mm of the sample is determined mechanically using a dial gauge.
As required from (39), three sequential measurements are taken at three different angles of incidence.In order to obtain a reliable determination of the permittivity data of the MUT, it is necessary to maximize the difference between the individual signal responses.This is the case when a measurement is made at the minimum and maximum possible angle, maximizing the difference in length of the signal path inside the MUT.The maximum usable angular range is however limited, but can be obtained from the considerations in Section III.For a 50.8 mm sized sample, the usable range is shown in Fig. 3. Therefore, transmission data at an incidence angle of 30 • , 45 • , and 60 • were taken.By adding further measurement data in the interval from 30 • to 60 • in 1 • -steps, as described in Section VI and VI-B, an improvement in the overall accuracy can be achieved.This is shown in Table I in terms of the RMSE to the mean value over the whole frequency range, since the alumina sample is expected to exhibit a nearlyconstant frequency response in the investigated frequency range.
The reference dielectric for the calibration procedure presented in Section VII, VII-A, and VII-B is a 0.2 mm thick slab of fused silica, with an assumed constant permittivity of 3.805 [37] in the frequency range under investigation.The angle of incidence for the calibration procedure Section VII and VII-A is chosen to be 45 • , and for Section VII and VII-B, the whole angular range of 30 • to 60 • is considered.The improvement that the inclusion of additional standards in Section VII, VII-A, and VII-B provide in comparison to Section VI, VI-A, and VI-B is shown in Table 1.Fig. 12 shows the extracted complex permittivity for both multivariate calibration methods.
To give a more reasonable comparison, especially since the absolute deviation of the complex permittivity depends strongly on the thickness of the MUT, the RMSE shall be expressed additionally in terms of S-parameters.Since the complex permittivity is determined directly using the method presented in Section VI, it must be converted into equivalent S-parameters for comparability with the other calibration methods.This is done by calculating (18) under normal incidence.
When comparing the RMSE of the S-parameter, the increased accuracy of the multivariate minimization methods compared with the univariate minimization methods is noticeable.However, this comes with the disadvantage of significantly higher number of measurements and thus measurement and computation time.The reason for the reduced RMSE with the minimization methods is essentially that besides averaging of noise, systematic error sources, which are not considered in the error model, are partially compensated by the inclusion of many measurements.This is especially true for phase drift of the VNA system, since it is differential for the two transmission measurements S 21 and S 12 and thus is partially compensated by the inclusion of both transmission directions.
Other systematic errors include a non-zero match and isolation error, uncertainties in the rotation angle and the determination of the Gaussian beam profile.These cause the loss tangent of the permittivity falling below zero, i.e. becoming unphysical, as can be seen in Fig. 12(b).However, the observed mean loss tangent of the reference measurement by the THz-TDS is within the uncertainty of the proposed calibration methods.It is expected, that a higher accuracy be in case of thicker material samples, as the contrast increases.

IX. CONCLUSION
In this work, a novel standardless calibration method for the characterization of dielectric materials has been presented, demonstrated by the characterization of a 0.625 mm thick alumina sample, measured in the frequency range from 220 GHz to 330 GHz.An extended calibration procedure using two known standards to characterize the same sample has also been presented.These calibration procedures require only a measurement of transmission for determination of the complex permittivity, demonstrating its usefulness for systems without combined transmit and receive modules.It has been shown that by measuring the MUT at an oblique angle of incidence and absorbing the signal reflected from the sample, the 12-term error matrix otherwise common to free-space setups is significantly simplified.This leaves only two unknown error terms in the forward case or three unknown error terms in the case including the backward direction.The preconditions necessary for this calibration and measurement, have been described in detail.Furthermore, a general description of beam propagation for collimated beams has been provided, which allows to accurately describe the problem of defocusing and lateral beam displacement present in typical highly-focused quasi-optical setups.The incorporation of higher-order Gaussian modes allows the use of horn antennas with low Gaussicity, enabling an easier and cheaper setup of the characterizing system.

FIGURE 1 .
FIGURE 1.(a) Schematic diagram of the symmetric free-space setup (top down view).The MUT (blue) is mounted rotatably.The center of rotation, the origin of the coordinate system and the focus point of the lenses are located on the surface of the MUT.(b) Measurement setup for the frequency range from 220 GHz to 330 GHz.

FIGURE 2 .
FIGURE 2. Signal flow graph of the simplified, transmission-only error model.

FIGURE 3 .
FIGURE 3. Measured estimation of isolation |E IF,est | and reflection match E M,est over the angle of incidence θ for the presented setup in Fig. 1.The usable angular range is marked yellow.

FIGURE 4 . 1 -
FIGURE 4. 1-port error model for the estimation of the reflection match.

FIGURE 5 .
FIGURE 5. Amplitude of the analytic transmission coefficient of 1 mm thick fused silica under different angles of incidence.The areas of A pos are marked yellow and A neg light yellow.

FIGURE 6 .
FIGURE 6. Magnitude profile of the transmission through a 1 mm thick fused silica sample (blue) and polynomial fit (red).

FIGURE 7 .
FIGURE 7. Optical path of an TE-polarized incident wave on a thin layer of a material surrounded by vacuum.

FIGURE 8 .
FIGURE 8. Shift and Fabry-Pérot echoes of the transmitted Gaussian beam E Tx inside a dielectric plate (red).The Gaussian beam E Rx associated with the receiver optics is depicted blue.Observable is the decreasing amplitude and beam broadening for consecutive reflected Gaussian beams.

FIGURE 9 .
FIGURE 9. Knife edge method.The metal sheet is angled to radiate the reflected signal out of the system.

FIGURE 10 .
FIGURE 10.Beam profile measurement at 275 GHz.(a) Measured (solid) and fitted (dashed) amplitude profile A(x).(b) Measured (solid) and fitted (dashed) Gaussian beam profile derived from the derivative of the amplitude profile of (a).

FIGURE 11 .
FIGURE 11.Gaussian beam parameters over the frequency.(a) Beam width.(b) Shift of beam waist.

FIGURE 12 .
FIGURE 12.Comparison of the extracted complex permittivity of 0.625 mm thick alumina.VI-B and VII-B refer to the methods proposed in their respective sections.(a) Real part of relative permittivity.(b) Loss tangent.