# Quasi-optical modeling of a millimeter- and submillimeter-wave free-space characterization bench

Gregory Gaudin, Clément Henry, Daniel Bourreau, Alain Peden

► To cite this version:

Gregory Gaudin, Clément Henry, Daniel Bourreau, Alain Peden. Quasi-optical modeling of a millimeter- and submillimeter-wave free-space characterization bench. International Journal of Microwave and Wireless Technologies, 2025, pp.1-11. 10.1017/S1759078724001314. hal-04958385v2

HAL Id: hal-04958385

https://imt-atlantique.hal.science/hal-04958385v2

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Quasi-Optical Modelling of a Millimetre- and Submillimetre-Wave Free-Space Characterisation Bench

Gregory Gaudin¹,²,³ and Clément Henry¹,²,⁴ and Daniel Bourreau¹,² and Alain Peden¹,²

¹ Microwaves department, IMT Atlantique, Brest, France

² Lab-STICC, UMR CNRS 6285, Brest, France

³ gregory.gaudin@imt-atlantique.fr

⁴ clement.henry@imt-atlantique.fr

## Abstract

In this article, we present a simulation tool for modelling a quasi-optical bench for material characterisation. The model uses a Gaussian beam expansion and tracking analysis, together with a modal analysis to enable comparison of the simulated reflection and transmission S-parameters to the ones measured with a 4-port Vector Network Analyser (VNA). A Thru-Reflect-Line (TRL) calibration is performed to de-embed the simulated S-parameters of a dielectric slab located between two lens antennas, showing good agreement with the analytical slab model used for experimental permittivity extraction.

## I Introduction

Quasi-optical (QO) techniques for characterising materials in free-space are particularly appealing as they allow for permittivity extraction over a wide frequency range, and are well-suited to millimetre and sub-millimetre wavelengths such as in the W (75-110 GHz) or J (220-330 GHz) bands [1, 2] and up to THz frequencies [3]. These techniques enable the characterisation of materials without specific machining of the sample, which is required when using rectangular or coaxial waveguide cells or cavities [4, 5]. To allow for precise permittivity extraction, the QO setup must be calibrated so that the actual measured S-parameters are compatible with the model employed in the extraction process (e.g. parallel-sided slab model under plane wave at normal incidence) [6]. The plane wave assumption is restricted to a local area where the beam and the material interact, as it is the case with a paraxial Gaussian beam, and may lead to an inaccurate model [7].

To evaluate the impact, on the extraction process, of any QO component misalignment during both the calibration and measurement steps or the impairments due to beams with non-planar wavefronts, modelling the test bench can be a valuable approach. Commercial software for full-wave electromagnetic simulations offer interesting capabilities but may be unsuitable for this specific context due to the large component sizes with respect to the wavelengths at these frequencies, leading to excessively long simulation times and a prohibitive memory cost. An intermediate solution involves using proprietary software like Feko,

which hybridizes full-wave and asymptotic methods, or GRASP which employs ray launching. However, these software do not provide trivially an end-to-end simulation pipeline for mapping the characterisation bench model to simulated S-parameters, and thus study the impact of multiple parameters such as misalignment and misplacement on the accuracy. Analytical models with Gaussian beam incidence can be also implemented to study the interaction with a slab and to extract its permittivity [7].

This paper addresses the modelling of a free-space characterisation bench with an in-house tool based on Gaussian beam expansion and tracking. Preliminary results of this work, without the complete modelling of the characterisation bench, have been presented at the French Microwave Days conference [8]. As performed during the experimental characterisation, a free-space Thru-Reflect-Line (TRL) calibration such as the one described in [9] is also applied to the simulated S-parameters.

The paper is structured as follows: Section II outlines the modelling approach. Section III provides a comparison between the simulation of a Gaussian optics lens antenna (GOLA) and near-field measurements obtained using a robotic arm. Section IV presents a simulation of the complete bench including a TRL calibration and a comparison of the results with a dielectric slab model. Finally, Section V covers the conclusion and future work.
## II Modelling of quasi-optical systems

### A Field expansion

When working with GOLAs, the Gaussian beam formalism [7, 10] is a convenient method for propagating the electromagnetic field. However, because the fields produced by some antennas are not purely Gaussian, it becomes practical to expand these fields into a set of Gaussian beams. There are several methods for decomposing a field into Gaussian beams. Instead of describing the radiated field only with the fundamental modes (TEM₀₁) of a Gaussian beam under the paraxial approximation, the field can be expanded into the Gauss-Laguerre modes (for circular aperture shapes), or Gauss-Hermite modes (for rectangular aperture shape) [10, 11]. However, all those modes share the same propagation axis, implying that the incoming field must exhibit weak divergence. A more suitable solution is to use an expansion surface to decompose the incoming field. A technique introduced in [12], known as multi-beam expansion, is used to decompose a field into a sum of vectorial Gaussian beams. This method efficiently expresses fields produced by a horn antenna feeding an optical system without requiring a high number of beams [13].

### 1 Vectorial Astigmatic Gaussian Beam

The modelling tool relies on the use of Vectorial [12] Astigmatic [14] Gaussian Beams. The electric field from such beams is expressed as

$$ ^{n}\mathbf{E}_{gb}(^{n}\boldsymbol{\rho}, z) = \begin{bmatrix} \boldsymbol{\alpha} \\ -\frac{1}{2}\boldsymbol{\alpha}^{\mathrm{T}}(\mathbf{Q}(z) + \mathbf{Q}^{\mathrm{T}}(z))\boldsymbol{\rho} \end{bmatrix} u(\boldsymbol{\rho}, z), \quad (1) $$

with $^{n}\mathbf{E}$ being the electric field written in frame (n), $^{n}\boldsymbol{\rho}$ the transverse coordinates written in frame (n), $\mathbf{Q} \in \mathbb{C}^{2 \times 2}$ the complex curvature matrix and $\boldsymbol{\alpha} \in \mathbb{C}^2$ a couple of transverse polarisation coefficients. The longitudinal component is derived from the partial derivatives of the transverse components [12].

In Eq. 1, $u$ is the scalar function [12, 14], defined as

$$ u(\boldsymbol{\rho}, z) = u_0(z) \exp \left\{ -j \frac{k}{2} \boldsymbol{\rho}^{\mathrm{T}} \mathbf{Q}(z) \boldsymbol{\rho} + j\phi(z) \right\}, \quad (2) $$

where $\boldsymbol{\rho}$ is the transverse position vector, and $z$ the longitudinal position, both written in the local reference frame (n) of the beam. The normalisation factor $u_0$ ensures that the beam carries 1 W through each cross-sectional area when propagating in a lossless medium.

The definitions of the complex curvature matrix $\mathbf{Q}$, the normalisation factor $u_0(z)$, and the complex phase $\phi(z)$ are provided in Table 1, where $k_0$ and $Z_0$ are the wave-number and the impedance in free-space, respectively.

### 2 Gaussian Beam Expansion

At each sampled point $\boldsymbol{r}_p$ on the expansion surface, the incoming electromagnetic field $\mathbf{E}_i(\boldsymbol{r}_p)$ is expanded into a sum of N Gaussian beams, originating from the N sampling points, as follows [12]:

$$ {}^{g}\mathbf{E}_{i}(\boldsymbol{r}_{p}) = \sum_{n=1}^{N} {}^{g}\mathbf{K}_{n}^{n} \mathbf{E}_{gb,n}(^{n}\boldsymbol{r}_{p}), \quad \forall p \in [1, N], \quad (4) $$

where ${}^{g}\mathbf{K}_n \in \mathbb{R}^{3 \times 3}$ represents the transformation matrix between the local frame (n) of the n-th beam and the reference frame (g) of the system. Each beam $\mathbf{E}_{gb,n}$ is divided into XZ and YZ polarisations,

$$ ^{n}\mathbf{E}_{gb,n} = \alpha_{x,n} \cdot {}^{n}\mathbf{E}_{xz,n} + \alpha_{y,n} \cdot {}^{n}\mathbf{E}_{yz,n}, \quad (5) $$

with $\alpha_{x,n}$ and $\alpha_{y,n}$ being the components of the polarisation vector $\boldsymbol{\alpha}$ of the n-th beam, mentioned in Eq. 1. Its propagation axis $\hat{z}$ is aligned with the Poynting vector of the incoming field at the corresponding sampled point $n$.

This leads to a linear system written in its matrix form

$$ \underbrace{\begin{bmatrix} \mathbf{E}_{xz \cdot x} & \mathbf{E}_{xz \cdot y} \\ \mathbf{E}_{yz \cdot x} & \mathbf{E}_{yz \cdot y} \end{bmatrix}}_{2N \times 2N} \underbrace{\begin{bmatrix} \mathbf{A}_x \\ \mathbf{A}_y \end{bmatrix}}_{2N \times 1} = \underbrace{\begin{bmatrix} \mathbf{E}_{s \cdot x} \\ \mathbf{E}_{s \cdot y} \end{bmatrix}}_{2N \times 1}, \quad (6) $$

where $\mathbf{A}_x$ and $\mathbf{A}_y$ are column vectors filled with the coefficients $\alpha_{x,n}$ and $\alpha_{y,n}$, respectively, and

$$ \mathbf{E}_{xz \cdot x} = \begin{bmatrix} \mathbf{E}_{xz,1}(\boldsymbol{r}_1) \cdot \hat{\mathbf{x}}_1 & \dots & \mathbf{E}_{xz,N}(\boldsymbol{r}_1) \cdot \hat{\mathbf{x}}_1 \\ \vdots & \ddots & \vdots \\ \mathbf{E}_{xz,1}(\boldsymbol{r}_N) \cdot \hat{\mathbf{x}}_N & \dots & \mathbf{E}_{xz,N}(\boldsymbol{r}_N) \cdot \hat{\mathbf{x}}_N \end{bmatrix}, \quad (7a) $$

and

$$ \mathbf{E}_{s \cdot x} = \begin{bmatrix} \mathbf{E}_s(\boldsymbol{r}_1) \cdot \hat{\mathbf{x}}_1 \\ \vdots \\ \mathbf{E}_s(\boldsymbol{r}_N) \cdot \hat{\mathbf{x}}_N \end{bmatrix}, \quad (7b) $$

for the XZ polarisations and local axes $\hat{\mathbf{x}}_n$. Note that similar definitions apply to the other terms of Eq. 6. Finally, solving 6 yields the vectors of N polarisation coefficients $\alpha_x$ and $\alpha_y$.

### B Gaussian beam tracking

The beams generated by the expansion are launched through the quasi-optical system. A tracking and bouncing operation is performed to model each beam interaction with the system components. It is assumed that for one incident beam, one reflected and one refracted beams are produced [12, 14, 15], each with its own polarisation coefficients and curvature matrix.

### 1 Computation of the polarisation coefficients

Assuming that the incoming beam spot size is small enough with respect to the surface at the intersection point, the polarisation coefficients of the reflected and
<caption>Table 1: Components of the Gaussian beam’s scalar function</caption>

$$
\mathbf{Q}(z) = \begin{bmatrix} \frac{\cos^2 \theta}{q_1 + z} + \frac{\sin^2 \theta}{q_2 + z} & \frac{1}{2} \sin 2\theta \left( \frac{1}{q_1 + z} + \frac{1}{q_2 + z} \right) \\ \frac{1}{2} \sin 2\theta \left( \frac{1}{q_1 + z} + \frac{1}{q_2 + z} \right) & \frac{\sin^2 \theta}{q_1 + z} + \frac{\cos^2 \theta}{q_2 + z} \end{bmatrix} \quad (3a)
$$

$$
\frac{Q_{12} + Q_{21}}{Q_{11} - Q_{22}} = \frac{\sin 2\theta}{\cos^2 \theta - \sin^2 \theta} = \frac{\sin 2\theta}{\cos 2\theta} = \tan 2\theta \quad (3b)
$$

$$
u_0(z) = \sqrt{\frac{k_0 Z_0}{\pi} \sqrt{4 \cdot \Im\{Q_{11}(z)\} \cdot \Im\{Q_{22}(z)\} - (\Im\{Q_{12}(z)\} + \Im\{Q_{21}(z)\})^2}} \quad (3c)
$$

$$
\phi(z) = -\gamma z + \frac{1}{2} \left( \angle \left\{ \frac{1}{q_1 + z} \right\} + \angle \left\{ \frac{1}{q_2 + z} \right\} \right) \quad (3d)
$$

refracted beam can be calculated using the Fresnel for-
mulas [14, 15]. The formulas need the polarisation co-
efficients to be expressed in a specific reference frame,
where each transverse axis represents either the (S) or
(P) polarisations. The S polarisation is perpendicu-
lar to the incident plane, while the P polarisation lies
within the incident plane. A transfer matrix, is used
to express the transverse polarisation $^i\alpha_I$ of the in-
coming beam in a new frame, referred to as (sp). The
characteristics of this frame, as well as the ones of re-
fracted (t) and reflected (r) beams, are detailed in [14].
The Fresnel formulas now apply to calculate the polar-
isation coefficients $^r\alpha_R$ of the reflected beam and $^t\alpha_T$
of the refracted beam

$$
{}^r \alpha_R = ({}^s \mathbf{K}_i {}^i \alpha_I) \odot \begin{bmatrix} \frac{n_i c_2 - n_t c_1}{n_i c_2 + n_t c_1} \\ \frac{n_i c_1 - n_t c_2}{n_i c_1 + n_t c_2} \end{bmatrix} \quad (8a)
$$

$$
{}^t \alpha_T = ({}^s \mathbf{K}_i {}^i \alpha_I) \odot \begin{bmatrix} \frac{2n_i c_1}{n_i c_2 + n_i c_1} \\ \frac{2n_i c_1}{n_t c_1 + n_i c_2} \end{bmatrix} \quad (8b)
$$

$$
c_1 = -\hat{n} \cdot \hat{z}_i, \quad c_2 = \sqrt{1 - \frac{n_i^2}{n_t^2}(1 - c_1^2)}, \quad (8c)
$$

where $\odot$ is the element-wise product, $\sqrt{a} = \sqrt{a}$ if
$a \ge 0$ and $\sqrt{a} = j\sqrt{|a|}$ otherwise, $\hat{n}$ the normal vec-
tor to the surface at the intersection point (oriented to-
ward the incoming medium), $n_i$ is the refractive index
of the incident medium, and $n_t$ is the refractive index
of the refracting medium.

## 2 Computation of the curvature matrices

A tangent frame to the surface, denoted as (uv), is
necessary to account for the curvature of the diop-
tre surface. The computation of the curvature matri-
ces of the generated beams is performed in this frame
before being expressed in the previously defined (r)
and (t) frames. The (uv) frame axes are $\hat{v} \equiv \hat{y}_i$ and
$\hat{u} = \hat{v} \times \hat{n}$.

To express the incident beam in the (uv) frame, ini-
tially defined in the (i) frame, we use transformation
matrices $^i\mathbf{K}_{uv}$ for the incident frame, $^r\mathbf{K}_{uv}$ for the re-
flected frame, and $^t\mathbf{K}_{uv}$ for the transmitted frame.

Next, the transfer matrices can be employed to define
the beam curvature matrices $^{uv}\mathbf{Q}_i$, $^{uv}\mathbf{Q}_r$, and $^{uv}\mathbf{Q}_t$
of the incident, reflected, and transmitted beams, respec-
tively. First, the curvature matrix $^a\mathbf{Q}_i$ (i.e. written in
frame (a)) reads

$$
^{uv} \mathbf{Q}_i = ({}^i \mathbf{K}_{uv})^{\mathrm{T}} {}^i \mathbf{Q}_i {}^i \mathbf{K}_{uv} \quad (9)
$$

The derivation of $^{uv}\mathbf{Q}_r$ and $^{uv}\mathbf{Q}_t$ can be found in [14]
and results in the following expressions

$$
^{uv} \mathbf{Q}_r = {}^{uv} \mathbf{Q}_i - 2 {}^{uv} \mathbf{C} (\hat{\mathbf{n}} \cdot \hat{\mathbf{z}}_i - \hat{\mathbf{n}} \cdot \hat{\mathbf{z}}_r) \quad (10a)
$$

$$
{}^t \mathbf{Q}_r = \left( {}^r \mathbf{K}_{uv} \right)^{-1} {}^{uv} \mathbf{Q}_r {}^r \mathbf{K}_{uv}^{-1} \quad (10b)
$$

$$
^{uv} \mathbf{Q}_t = \frac{n_i}{n_t} \left\{ ^{uv} \mathbf{Q}_i - {}^{uv} \mathbf{C} \left( \hat{\mathbf{n}} \cdot \hat{\mathbf{z}}_i - \frac{n_t}{n_i} \hat{\mathbf{n}} \cdot \hat{\mathbf{z}}_t \right) \right\} \quad (10c)
$$

$$
{}^t \mathbf{Q}_t = ({}^t \mathbf{K}_{uv})^{-1} {}^{uv} \mathbf{Q}_t {}^t \mathbf{K}_{uv}^{-1} \quad (10d)
$$

where $^{uv}\mathbf{C} \in \mathbb{R}^{2 \times 2}$ is the surface curvature matrix at the intersection point.

When the surface is curved, the generated beams can
no longer be simple astigmatic [14, p.6031]. A gen-
eral method to retrieve the parameters $q_1$, $q_2$, and $\theta$
required for beam propagation (Eq. 3a) is needed to
consider them as generally astigmatic. The parameters
$q_1$ and $q_2$ will be the inverses of the eigenvalues of $\mathbf{Q}$
[14, p.6032] and the complex angle $\theta$ can be calculated
using Eq. 3b.

C Matching of the generated beams

In the studies conducted in [14, 15], the Gouy phase
$\eta$ and the accumulated phase $\phi_{ac}$ are neglected in
the calculations. However, it is essential to account
for them to avoid introducing phase discontinuities.
Moreover, the amplitude of each beam is described in
relation to the position of their respective waists. Since these positions are not aligned, this also requires an adjustment in amplitude. Therefore, a coupling term $a_0 \exp\{j\phi_0\}$ must be added:

$$a_{0,g} = \frac{|u_i(\mathbf{0}, d_i)|}{u_{0,g}(\mathbf{0}, 0)}, \quad g \in \{r, t\}, \quad (11a)$$

$$\phi_{0,g} = \angle u_i(\mathbf{0}, d_i) - \eta_g(0), \quad g \in \{r, t\}. \quad (11b)$$

where $d_i$ is the distance between the intersection and the incident beam origin. The prefix $i$ identifies the incident beam ($i$), while $g$ is for either refracted ($t$) or reflected ($r$).

### D Recombination and probe modelling

The electrical field at a point $r$ is obtained using the recombination formula given by Eq. 12

$$^{g}\mathbf{E}(^{g}\mathbf{r}) \approx \sum_{n=1}^{N} {}^{g}\mathbf{K}_{n}^{n} \mathbf{E}_{gb,n}(^{n}\mathbf{r}) \quad (12)$$

To account for the presence of the probe, the received fields $E$ and $H$ at its aperture $A_p$ are expressed as a combination of guided modes [16, p.120],

$$\mathbf{E} = \sum_{N} c_k e_k. \quad (13)$$

The modal coefficients $c_k$ are derived using the orthonormal properties of these modes. This yields

$$c_k = \iint_{A_p} \mathbf{E} \times \mathbf{h}_k^* \cdot \hat{z} dS \quad (14)$$

where $e_k$ and $h_k$ denote the electric and magnetic fields of the guided mode $k$. A mode-matching technique provides a generalised scattering matrix $S_{gsm}$ for computing the modal coefficients received by the waveguide feeding the probe [16, p. 122], such as

$$\begin{bmatrix} \sim \\ \mathbf{b}_2 \end{bmatrix} = \mathbf{S}_{GSM} \begin{bmatrix} c \\ \mathbf{0}_N \end{bmatrix}, \quad (15)$$

where $\mathbf{0}_N$ is a column vector of N zeros. These coefficients allow for the determination of transmission parameters $b_2$ for each mode, giving $S_{21} = \text{diag}(b_2)$, considering that the source emits a total power of 1 W. In the case of an Open Ended Waveguide (OEWG) probe, this matrix can be neglected, as it only induces a phase shift. Employing the same methodology with a probe located at the source and considering only backward-propagating waves enables retrieving the $S_{11}$ parameter.

## III Application to a GOLA

### A Comparison with measurements done with a robotic arm

The objective of this section is to validate experimentally the results obtained using the simulation tool. An

experimental setup in the J-band was implemented to carry out the characterisation of the electromagnetic field radiated by a diagonal horn associated with a lens.

### 1 Measurement setup

The experimental setup is divided into two parts. The first part remains fixed and includes a millimetre-wave transmitting/receiving converter connected to a diagonal horn (referred to as the source) and a PTFE lens. The second part consists of a UR10e robotic arm with six degrees of freedom, designed by Universal Robots (visible in Fig. 1). This arm is controlled by its digital counterpart through the RobotDK software and receives instructions via a Python script. A second transmitting/receiving converter with an OEWG probe is installed on the robotic arm. A Vector Network Analyser (VNA) is used to measure the transmission parameter ($S_{21}$) of the entire system. Before acquiring data, the setup is calibrated using the well-known Thru-Reflect-Line guide-to-guide method [17], the horns on the converters being removed. Although it is common practice to apply corrections to measurements obtained from an OEWG probe, no additional post-processing has been applied since the measured data are directly compared to a simulation that accounts for the presence of the probe.

### 2 Measurements within the E-Plane

The probe scans the YZ-plane, first along the y-axis and then towards the z-axis. This approach mitigates the occurrence of ripples, as observed in Fig. 2, where the measurements were performed with low spatial resolution and with no lens in front of the horn. This phenomenon can be explained by the rapid phase variation along the propagation axis, compared to the slower variation in the transverse dimensions. Consequently, any misalignment of the probe will have a more detrimental impact if the scan is performed first along the most sensitive dimension. Fig. 3 illustrates the results of a scan on the YZ plane with a step size of $1\lambda$ on the y-axis and $0.1\lambda$ on the z-axis and with a lens

<figure><img src="image_2.png" /><figcaption>Figure 1: Measurement setup for probing the field emitted by a lens antenna in the J-band, using a robotic arm.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Impact on the phase measured by the probe at 300 GHz with respect to the robotic arm movement scheme.</figcaption></figure>

placed at a 90 mm distance from the horn antenna. The measurements and the simulation are in good agreement in terms of magnitude and phase, demonstrating an almost perfectly paraxial output beam as required in the material characterisation setup.

## IV Modelling of the characterisation bench

The aim of this section is to demonstrate the capability of the simulation tool to perform a simulated free-space Thru-Reflect-Line (TRL) calibration on the S-parameters retrieved from simulation, to mimic the procedure conducted on the characterisation test bench, which consists of a material under test almost centrally positioned between two GOLAs.

### A Model setup

#### 1 Simulation

The previously described and implemented tool is used for modelling the setup depicted in Fig. 4, where GBE stands for Gaussian Beam Expansion. As a beam propagates, it diverges and covers a larger area when intersecting a surface. When the setup elements are not highly coupled, accuracy is maintained by performing a beam summation and field reexpansion right before the element to be modelled. In this simulation setup, four expansion surfaces are established. Ports 1 and 2 are rectangular apertures at which a modal analysis is performed to retrieve the reflection and transmission parameters $S_{11}$ and $S_{21}$, respectively. The radiation from port 1 is modelled as a Gaussian beam with a 2.11 mm waist. Lenses 1 and 2 are placed at their front focal lengths from their respective ports (i.e. horn aperture). The lenses and

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison between measurement (left) and simulation (right), of a GOLA operating at 300 GHz, as measured with a probe on the YZ plane.</figcaption></figure>

the slab are obtained from CAD models imported into the simulator with Gmsh [18], allowing the creation of meshes suitable for beam tracing.

Fig. 5 depicts the pipeline used to model the characterisation bench, considering the following assumptions:

* The reflection of the beam on lens 2 does not contribute to the reflection coefficient $S_{11}$.

* The reflection of the backward beam on lens 1 does not contribute to the transmission coefficient $S_{21}$.

* The lenses are lossless (zero extinction coefficient). Since this loss is present in both calibration and measurement, it is assumed to have no significant impact.

* The power received by the fundamental mode at the aperture is fully transferred into the waveguide.

* Both lenses are perfectly identical in terms of shape and relative permittivity.

#### 2 Slab analytical model

An infinite dielectric slab made of rexolite with parallel faces can be modeled using transmission line theory. Given $\epsilon_r$ as the permittivity of the slab, the reflection coefficient of a single interface air-rexolite is

$$ \Gamma = \frac{1 - \sqrt{\epsilon_r}}{1 + \sqrt{\epsilon_r}}, \qquad (16) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Simulation setup.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Simulation pipeline, where Σ stands for beam summation and G.B.E. for Gaussian Beam Expansion.</figcaption></figure>

and with γ as the propagation constant in the rezo-lite, the reflection and transmission coefficients of the dielectric slab under normal incidence can be determined:

$$S_{11,th} = \frac{\Gamma(1 - 2 \exp \{-2j\gamma h\})}{1 - \rho^2 \exp \{-2j\gamma h\}}, \quad (17a)$$

$$S_{21,th} = \frac{(1 - \Gamma^2) \exp \{-2j\gamma h\}}{1 - \rho^2 \exp \{-2j\gamma h\}}. \quad (17b)$$

For a slab made of an isotropic material, reciprocity is applied, meaning $S_{22,th} = S_{11,th}$ and $S_{12,th} = S_{21,th}$.

### 3 Thru-Reflect-Line calibration

A free-space TRL calibration is carried out on the characterisation bench for de-embedding the S-parameters measured with the VNA and thus extract the ones of the slab, a crucial step for accurately extracting its permittivity [6]. The procedure involves three steps, denoted T, R, and L and is performed only for the direction Port 1 → Port 2, under the assumption of perfect symmetry in the simulated system. Therefore, it is assumed that all optical elements are perfectly aligned, representing the ideal case for the bench. We propose here to perform this calibration on the simulated data obtained from the whole setup and compare the results with the dielectric slab theoretical model. In the present study, the element 'OBJ' shown in Fig. 4 is substituted by a perfect reflector for step R, an air line of about $λ/4$ (where λ is the mid-band wavelength) for step L, and then no object for step T. This process is synthesised in Fig. 6.

During a slab measurement, the GOLA at port 2 is not moved and remains in the configuration used for the

Thru step. Therefore, in order to compare the simulation to the analytical model of the slab, $S_{22,th}$, $S_{12,th}$ and $S_{21,th}$ must be adjusted to the correct reference plane (i.e., the face toward the port 1 of the slab here), giving

$$S_{22,th,c} = S_{22,th} \cdot \exp \{2jk_0h\}, \quad (18a)$$

$$S_{21,th,c} = S_{21,th} \cdot \exp \{jk_0h\}. \quad (18b)$$

To eliminate any offset caused by misplacement of either the slab or one of the two GOLAs with respect to the reference plane of the TRL, the phase of the product between $S_{11}$ and $S_{22}$ is utilised to compare with the model.

For each step, results are stored in S2P-format files to be processed by a MATLAB code enabling a TRL calibration.

## B Results

### 1 Simulation of the bench

The simulation involved a frequency sweep from 220 GHz to 330 GHz (over 10 points per GHz) for the three steps (T, R, and L depicted in Fig. 7), followed by the simulation with the slab (defined as step M). The results for the M step are presented in Fig. 8, which shows the simulated S-parameters obtained with and without considering internal reflections within the lenses, as well as the corrected analytical slab S-parameters. Resonances are accurately identified, and the S-parameter magnitudes from both sets exhibit similar trends. The ripples observed in Fig. 8a are attributable to multi-reflection inside the lenses, as evidenced by the smoothness of curve in Fig. 8b, which represents the simulation conducted without internal reflections in the lenses.

For the frequency sweep involving the four steps (T, R, L, and M), the process typically lasts about 8 hours on

<figure><img src="image_7.png" /><figcaption>Figure 6: Principle of the Thru-Reflect-Line calibration applied to the characterisation bench for de-embedding the S-parameters.</figcaption></figure>
a 12-core Intel Core i7-12700H CPU, and requires less
than 10 GB of RAM. A frequency sweep involving
only the step M spans less than 2 hours.

<figure><img src="image_8.png" /><figcaption>Figure 7: Magnitude of simulated $S_{11}$ (red) and $S_{21}$ (blue) of the bench at 220-330 GHz for the Thru-Reflect-Line steps of the calibration.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(a) Simulation with multi-reflections inside lenses.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Simulation without multi-reflections inside lenses.</figcaption></figure>

<caption>Figure 8: Magnitude of raw simulated S-parameters of the bench and analytical S-parameters of the slab at 220-330 GHz with a 12.815 mm rexolite slab ($\epsilon_r = 2.5329 - 0.00667j$). It is assumed that $S_{22} = S_{11}$ and $S_{12} = S_{21}$.</caption>

## 2 Application of the TRL

A TRL calibration is applied to the simulated S-
parameters and presented in Fig. 9. The figure demon-
strates a strong agreement between the extracted S-
parameters ($S_{22,s} \cdot S_{11,s}$ and $S_{21,s}$), and the cor-
rected analytical S-parameters ($S_{22,th,c} \cdot S_{11,th,c}$ and
$S_{21,th,c}$). Additionally, there is a strong correlation
with the measurements of $S_{11,m}$ and $S_{21,m}$ obtained
from the characterisation bench. This suggests that us-
ing a large beam or a plane wave yields nearly equiva-
lent results under normal incidence.

<figure><img src="image_11.png" /><figcaption>Figure 9: De-embedded reflection and transmission parameters of the dielectric slab obtained from the simulation, compared to the analytical model and to the measurements, with $\epsilon_r = 2.5329 - 0.00667j$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: De-embedded reflection and transmission parameters of the slab ($\epsilon_r = 2.5329 - 0.00667j$). GOLA lenses were shifted horizontally by 15 mm, the slab positioned 150 mm away from the beam waist, and the distance between ports 1 and 2 increased to 1 m.</figcaption></figure>

Next, we study the impact of a misplacement in the
characterisation bench. Figure 10 shows simulations
and measurements with a 15 mm lens horizontal off-
set from the horn apertures, causing a slightly more
focused beam. Although the slab is positioned away
from the beam waist, simulations and measurements
still well agree, and extracted permittivity remained
accurate.

## V CONCLUSION

A numerical tool has been implemented to simulate
quasi-optical systems, such as a free-space characteri-
sation bench, at millimetre and sub-millimetre wave-
lengths. This tool has been employed to simulate
a GOLA, a key component of the bench, showing
good agreement with measurements obtained using a
robotic test bench.

Additionally, it has been demonstrated that this tool
can effectively model the entire characterisation bench
and perform the slab S-parameter de-embedding nec-
essary for permittivity extraction. Future works in-
volve simulating a more realistic bench, accounting for
setup imperfections like tilts and misalignments to as-
sess their impact on permittivity extraction accuracy.

Acknowledgment

The work received support from France 2030 in
the PEPR (Priority Research and Equipment Pro-
gramme) programme "Future networks", and the
NF-SYSTERA project funded by the French ANR
(Agence Nationale de la Recherche) under Grant 22-
PEFT-0006.

Competing interests

The author(s) declare none.

Bibliography

[1] Óscar García Pérez, F. Tercero, and Samuel Lopez Ruiz. Free-space W-band setup for the electrical characterization of materials and mm-wave components. Technical report, Observatorio de Yebes Guadalajara (Spain), 06 2017.

[2] Elena Saenz, Luis Rolo, Kees Van't Klooster, Maurice Paquay, and Vladimir V. Parshin. Accuracy assessment of material measurements with a quasi-optical free-space test bench. In *6th European Conference on Antennas and Propagation (EuCAP)*, pages 572–576, 03 2012. ISBN 978-1-4577-0918-0. doi: 10.1109/EuCAP.2012.6206324.

[3] Jonathan Hammler, Andrew Gallant, and Claudio Balocco. Free-space permittivity measurement at terahertz frequencies with a vector network analyzer. *IEEE Transactions on Terahertz Science and Technology*, **6**, 10 2016. doi: 10.1109/TTHZ.2016.2609204.

[4] Hong Eun Choi, Wonjin Choi, Evgenya I. Simakov, Muhammed Zuboraj, Bruce E. Carlsten, and EunMi Choi. Error tolerant method of dielectric permittivity determination using a te01 mode in a circular waveguide at the W -band. *IEEE Transactions on Microwave Theory and Techniques*, **68** (2):808–815, 2020. doi: 10.1109/TMTT.2019.2951156.

[5] Nicolas Pantano, John Slabbekoorn, Fabrice Duval, and Eric Beyne. Broadband permittivity characterization of polymers up to 110 ghz using co-planar waveguides. In *2021 IEEE 71st Electronic Components and Technology Conference (ECTC)*, pages 1252–1257, 2021. doi: 10.1109/ECTC32696.2021.00203.

[6] Daniel Bourreau, Alain Peden, and Sandrick Maguer. A quasi-optical free-space measurement setup without time-domain

gating for material characterization in the W-band. *Instrumentation and Measurement, IEEE Transactions on*, 55:2022-2028, 01 2007. doi: 10.1109/TIM.2006.884283.

[7] Michal Mrnka, Roger Appleby, and Elena Saenz. Accurate S-parameter modeling and material characterization in quasi-optical systems. *IEEE Transactions on Terahertz Science and Technology*, 12(2):199-210, 2022. doi: 10.1109/TTHZ.2021.3140201.

[8] Gregory Gaudin, Daniel Bourreau, Clément Henry, and Alain Peden. Modélisation de systèmes quasi-optiques dédiés à la caractérisation de matériaux et à l'imagerie en bande J (220-330 GHz). In *JNM 2024, Antibes Juan-Les-Pins, France*, June 2024. URL https://imt-atlantique.hal.science/hal-04631622.

[9] Daniel Bourreau and Alain Péden. Solid and non-solid dielectric material characterization for millimeter and sub-millimeter wave applications. In *2020 50th European Microwave Conference (EuMC)*, pages 909-912, 2021. doi: 10.23919/EuMC48046.2021.9338045.

[10] Paul F. Goldsmith. Gaussian Beam Propagation, pages 9–38. Wiley-IEEE Press, 1998. ISBN 978-0-470-54629-1. doi: 10.1109/9780470546291.ch2.

[11] *Beam Optics*, chapter 3, pages 80–107. John Wiley Sons, Ltd, 1991. ISBN 9780471213741. doi: https://doi.org/10.1002/0471213748.ch3. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/0471213748.ch3.

[12] Alexandre Chabory, Jérôme Sokoloff, Sylvain Bolioli, and Paul François Combes. Computation of electromagnetic scattering by multilayer dielectric objects using gaussian beam based techniques. 6(6):654–662. ISSN 16310705. doi: 10.1016/j.crhy.2005.06.011.

[13] Alexandre Chabory, Jerome Sokoloff, and Sylvain Bolioli. Physically based expansion on conformal gaussian beams for the radiation of curved aperture in dimension 2. *Microwaves, Antennas Propagation, IET*, 2:152 – 157, 04 2008. doi: 10.1049/iet-map:20060168.

[14] Evgenia Kochkina, Gudrun Wanner, Dennis Schmelzer, Michael Tröbs, and Gerhard Heinzel. Modeling of the general astigmatic gaussian beam and its propagation through 3D optical systems. *Applied Optics*, 52(24):6030–6040, Aug 2013. doi: 10.1364/AO.52.006030.

[15] Arash Rohani, Amir Ahmad Shishegar, and S. Safavi-Naeini. A fast gaussian beam tracing method for reflection and refraction of general vectorial astigmatic gaussian beams from general curved surfaces. *Optics Communications*, 232:1–10, 03 2004. doi: 10.1016/j.optcom.2003.11.044.

[16] Jorge Ruiz-Cruz, José Montejo-Garai, and Jesus Rebollar. *Computer Aided Design of Waveguide Devices by Mode-Matching Methods*. InTech, 04 2010. ISBN 978-953-307-083-4. doi: 10.5772/9403.

[17] Glenn F. Engen and Cletus A. Hoer. Thru-reflect-line: An improved technique for calibrating the dual six-port automatic network analyzer. *IEEE Transactions on Microwave Theory and Techniques*, 27(12):987-993, 1979. doi: 10.1109/TMTT.1979.1129778.

[18] Christophe Geuzaine and Jean-François Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79(11):1309–1331, September 2009. ISSN 00295981. doi: 10.1002/nme.2579.