*Article* **Numerical Analysis of the Influence of Fabrication Process Uncertainty on Terahertz Metasurface Quality**

**Przemyslaw Lopato 1,\* , Michal Herbko 1,\* , Paulina Gora <sup>1</sup> , Ulrich Mescheder <sup>2</sup> , Andras Kovacs <sup>2</sup> and Alexander Filbert <sup>2</sup>**


**Abstract:** The purpose of this paper was to investigate the influence of fabrication process uncertainty on terahertz metasurface quality. The focus was on the effect of metasurface fabrication inaccuracy on resonances. To the best of our knowledge, this is the first paper to study the effect of the metasurface fabrication process on its resonant frequency. The terahertz split ring resonator-based metasurface is under consideration. Using a numerical model, the influence of the uncertainty of various geometrical parameters obtained during the fabrication process (mainly layer deposition, photolithography, and etching processes) is analyzed according to the resonance of the designed metasurface. The influence of the following parameters causes a shift of resonant frequencies of the considered metasurface: etching deviation *e*, metallization thickness *t*Al and SiO<sup>2</sup> layer thickness *t*SiO2. The quality of the metasurface affected by the variations of obtained geometrical parameters was determined by the deviation of resonant frequency ∆*f*r. The developed numerical model was verified by THz-TDS (terahertz time-domain spectroscopy) measurements of the fabricated structure.

**Keywords:** metasurfaces; terahertz devices; electromagnetic waves; finite element model

## **1. Introduction**

A metamaterial (MM) is an artificial material/device consisting of an array of subwavelength-size structural elements (shown in Figure 1) that can take on the properties of non-naturally occurring materials [1]. It can be more specifically described as a *ε*-negative (ENG) or *µ*-negative (MNG) material of either negative permittivity or permeability, respectively, a double-negative (DNG) when both of these quantities are simultaneously negative in a particular frequency range, and *ε*-near-zero (ENZ) or a *µ*-near-zero (MNZ) material with small, close-to-zero values for these quantities (this might be the case when the materials are close to their electric or magnetic plasma frequency) [2]. This creates new opportunities to interact with electromagnetic waves (from microwave band to visible light) and mechanical (ultrasonic) waves. In this paper, we will focus on electromagnetic waves in the terahertz frequency band. In this case, various field transformations of arbitrary incident fields can be achieved. Planar versions of MM-metasurfaces (MS) are much easier to fabricate (in comparison to 3D structures) and can potentially have many applications. The first of the areas in which metasurfaces can be used is telecommunications. Among other things, it can be used to improve the efficiency of antennas by increasing directivity, increasing power gain, increasing selectivity, or by multiplexing signals [3–6]. Other peripheral areas for metasurfaces are selective absorbing structures [7,8], polarizers [9,10], surface wave control devices [11,12], wave front manipulation devices [13,14] and imaging systems (for lensing electromagnetic waves) [15,16]. Particularly noteworthy are works showing plasmonic metasurfaces used as perfect absorbers and pulse controllers [17,18].

**Citation:** Lopato, P.; Herbko, M.; Gora, P.; Mescheder, U.; Kovacs, A.; Filbert, A. Numerical Analysis of the Influence of Fabrication Process Uncertainty on Terahertz Metasurface Quality. *Electronics* **2023**, *12*, 2198. https://doi.org/10.3390/ electronics12102198

Academic Editors: Raed A. Abd-Alhameed, Naser Ojaroudi Parchin, Mohammad Ojaroudi and Giovanni Crupi

Received: 2 March 2023 Revised: 5 May 2023 Accepted: 10 May 2023 Published: 12 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Other peripheral areas for metasurfaces are measurement systems. Among the most frequently mentioned are applications for evaluating biological samples, where it is most often indicated in the detection of biomolecules, microorganisms, cancer, or tumors [19]. Other applications include testing chemicals and measuring strain [20,21]. Exemplary metasurface structures are presented in Figure 1. signed structures resulting from the inaccuracy of microfabrication. This problem becomes particularly important in the higher bands of the spectrum—infrared and visible light—where (as the electromagnetic wavelength decreases) the inaccuracies in the manufacturing process increase significantly. The developed numerical model was verified by THz-TDS measurements of the fabricated structure.

*Electronics* **2023**, *12*, x FOR PEER REVIEW 2 of 13

plary metasurface structures are presented in Figure 1.

frequency.

imaging systems (for lensing electromagnetic waves) [15,16]. Particularly noteworthy are works showing plasmonic metasurfaces used as perfect absorbers and pulse controllers [17,18]. Other peripheral areas for metasurfaces are measurement systems. Among the most frequently mentioned are applications for evaluating biological samples, where it is most often indicated in the detection of biomolecules, microorganisms, cancer, or tumors [19]. Other applications include testing chemicals and measuring strain [20,21]. Exem-

Characterization of the fabrication process uncertainty is a key question for quality control. So far, the effect of changing the geometry of the metasurface created by inaccuracies in the fabrication process has been studied for different dimensions of rectangular structural elements in the work [22]. In our work, we investigated the effect of the fabrication process on the operation of SRR-based meta-surfaces; it is a more complicated structure to analyze because as a result of the process, some dimensions increase (gap width and distances between elements) and some decrease (line width and size of the square loop) as a result of over-etching. In our work, the effect of not only changing the geometry as a result of the etching process but also of oxidation and metallization was studied. The paper [22] studied the effect of the fabrication process on the metasurface in a hologram application, where the effect of fabrication on the phase of the electromagnetic wave was examined. In this paper, the effect of the fabrication process on resonances was studied, as we are interested in sensor applications. To the best of our knowledge, this is the first paper to study the effect of the fabrication process of a metasurface on its resonant

The fabrication of the MS in the terahertz range is carried out using a silicon-based microtechnology process. Although this technology is precise, there are some deviations from the planned geometric dimensions or thicknesses of individual layers. In this paper, the split-ring Resonator (SRR)-based metasurface was designed for the resonant frequency of *f*<sup>r</sup> = 1 THz, and the influence of fabrication process uncertainty on terahertz metasurface quality was characterized using numerical analysis. In this article, the square version of the SRR was chosen. The quality of the metasurface affected by the received variation geometrical parameters was determined by the resonant frequency shift Δ*f*r. The increase in the metallization thickness *t*Al as well as the SiO2 layer thickness *t*SiO2 causes an increase in the resonant frequencies of the considered metasurface, while etching deviation *e* has the opposite effect. The ±20% parameter variation range was assumed. The resulting frequency deviations potentially causing a degradation of metasurface performance in some applications were examined. Thanks to the use of a validated numerical model, it was possible to estimate the deviations of electromagnetic parameters of the de-

**Figure 1. Figure 1.**Exemplary metasurface structures. Exemplary metasurface structures.

Characterization of the fabrication process uncertainty is a key question for quality control. So far, the effect of changing the geometry of the metasurface created by inaccuracies in the fabrication process has been studied for different dimensions of rectangular structural elements in the work [22]. In our work, we investigated the effect of the fabrication process on the operation of SRR-based meta-surfaces; it is a more complicated structure to analyze because as a result of the process, some dimensions increase (gap width and distances between elements) and some decrease (line width and size of the square loop) as a result of over-etching. In our work, the effect of not only changing the geometry as a result of the etching process but also of oxidation and metallization was studied. The paper [22] studied the effect of the fabrication process on the metasurface in a hologram application, where the effect of fabrication on the phase of the electromagnetic wave was examined. In this paper, the effect of the fabrication process on resonances was studied, as we are interested in sensor applications. To the best of our knowledge, this is the first paper to study the effect of the fabrication process of a metasurface on its resonant frequency.

The fabrication of the MS in the terahertz range is carried out using a silicon-based microtechnology process. Although this technology is precise, there are some deviations from the planned geometric dimensions or thicknesses of individual layers. In this paper, the split-ring Resonator (SRR)-based metasurface was designed for the resonant frequency of *f*<sup>r</sup> = 1 THz, and the influence of fabrication process uncertainty on terahertz metasurface quality was characterized using numerical analysis. In this article, the square version of the SRR was chosen. The quality of the metasurface affected by the received variation geometrical parameters was determined by the resonant frequency shift ∆*f*r. The increase in the metallization thickness *t*Al as well as the SiO<sup>2</sup> layer thickness *t*SiO2 causes an increase in the resonant frequencies of the considered metasurface, while etching deviation *e* has the opposite effect. The ±20% parameter variation range was assumed. The resulting frequency deviations potentially causing a degradation of metasurface performance in some applications were examined. Thanks to the use of a validated numerical model, it was possible to estimate the deviations of electromagnetic parameters of the designed structures resulting from the inaccuracy of microfabrication. This problem becomes particularly important in the higher bands of the spectrum—infrared and visible light—where (as the electromagnetic wavelength decreases) the inaccuracies in the manufacturing process increase significantly. The developed numerical model was verified by THz-TDS measurements of the fabricated structure.

#### **2. THz Metasurface Design and Fabrication**

The arbitrary selected split-ring resonator metasurface was designed with the following assumptions: a primary resonant frequency of *f*r1 = 1 THz, and a secondary resonant frequency of *f*r2 < 3 THz. Moreover, we assumed that the magnitude of transmission

coefficient obtained within simulations |*S*21| < −40 dB for both resonances. Additionally, there were some restrictions according to the planned fabrication process that will be described later in this paper. The design (determination of dimensions) was based on a simple optimization procedure. In the first step based on the literature review, an initial set of geometrical parameters was developed. Relationships from the paper [23] were used to develop the initial design. Then, this model was implemented and simulated in the COM-SOL Multiphysics environment. In the third step, the geometrical parameters describing the structural element were tuned to meet the assumed electromagnetic parameters (*f*r1, *f*r2, and |*S*21|). The obtained dimensions are presented in Table 1 and the geometry of the squared SRR structural element with the description of dimensions is shown in Figure 2. The structure was developed on a 520 µm thick silicon (Si) wafer. First, the wafer was covered with a thermal oxide layer (SiO2) (on both sides) with a thickness of 170 nm. The front side of the wafer was coated with thermally evaporated aluminum. The wafer was coated with an AZ 1518 photoresist layer and exposed with the designed terahertz filter mask. The aluminum layer was wet-chemically structured with a phosphoric acid etching mixture, and the photoresist layer was removed in an oxygen plasma. The individual chips for the measurement were separated with a diamond saw. The flowchart of the manufacturing process is shown in Figure 3. The obtained line width was 2.5 µm, while the capacitive gap was 2.5 µm and the length of the side of the square loop was 21.7 µm. *Electronics* **2023**, *12*, x FOR PEER REVIEW 3 of 13 **2. THz Metasurface Design and Fabrication** The arbitrary selected split-ring resonator metasurface was designed with the following assumptions: a primary resonant frequency of *f*r1 = 1 THz, and a secondary resonant frequency of *f*r2 < 3 THz. Moreover, we assumed that the magnitude of transmission coefficient obtained within simulations |*S*21| < −40 dB for both resonances. Additionally, there were some restrictions according to the planned fabrication process that will be described later in this paper. The design (determination of dimensions) was based on a simple optimization procedure. In the first step based on the literature review, an initial set of geometrical parameters was developed. Relationships from the paper [23] were used to develop the initial design. Then, this model was implemented and simulated in the COM-SOL Multiphysics environment. In the third step, the geometrical parameters describing the structural element were tuned to meet the assumed electromagnetic parameters (*f*r1, *f*r2, and |*S*21|). The obtained dimensions are presented in Table 1 and the geometry of the

> **Table 1.** Dimensions of designed SRR metasurface. The structure was developed on a 520 µm thick silicon (Si) wafer. First, the wafer was


squared SRR structural element with the description of dimensions is shown in Figure 2.

**Figure 2.** Dimensions of SRR metasurface on the silicon wafer with a top view (**top**) and cross-section (**bottom**). **Figure 2.** Dimensions of SRR metasurface on the silicon wafer with a top view (**top**) and crosssection (**bottom**).

**Figure 3.** Flowchart of the manufacturing process. **Figure 3.** Flowchart of the manufacturing process.

**Table 1.** Dimensions of designed SRR metasurface. **Dimensions [µm]** Periodicity, *p* 24.2 Square loop size, *a* 21.7 Line width, *w* 2.5 Gap, *g* 2.5 Aluminum thickness, *t*Al 0.5 The design and further analyses were carried out using a finite element method (FEM) based numerical model. All the computations were made in the COMSOL Multiphysics 6.0 software. The geometry of the model is shown in Figure 4. The following values of dielectric properties were assumed for the materials. The refractive index of silicon is *n*Si = 3.48 (permittivity is *ε*rSi = 12.11) [24]. The real and imaginary parts of the complex permittivity of silicon dioxide for the terahertz range are *ε* 0 rSiO2 = 4.45 and *ε* 00 rSiO2 = 0.03, respectively [25]. The Drude model was used to describe the value of the aluminum refractive index *nAl* [26]: in the case of frequency domain computations, some sophisticated methods for its reduction should be applied. In this high-frequency electromagnetic problem the following frequency domain equation is solved [28,29]: ∇ൈ୰ ିଵሺ∇ൈሻ െ <sup>ଶ</sup> ൬୰ െ ൰ ൌ 0 (2) where: *E*—electric field, *µ*r—relative permeability, *ε*0—permittivity of vacuum, *ε*r—rela-

frequency range (multiple, relatively high amplitude ripples obtained in spectrum); thus,

$$m\_{Al} = \sqrt{\varepsilon\_{rAl}^{\prime} + i\varepsilon\_{rAl}^{\prime}} = \sqrt{1 - \frac{\omega\_p^2}{\omega^2 + \gamma^2} + i\frac{\omega\_p^2 \gamma}{\omega\_p^3 + \gamma^2 \omega}}\tag{1}$$

(1)

(FEM) based numerical model. All the computations were made in the COMSOL Multiphysics 6.0 software. The geometry of the model is shown in Figure 4. The following values of dielectric properties were assumed for the materials. The refractive index of silicon is *n*Si = 3.48 (permittivity is *ε*rSi = 12.11) [24]. The real and imaginary parts of the complex permittivity of silicon dioxide for the terahertz range are *ε*′rSiO2 = 4.45 and *ε*″rSiO2 = 0.03, where *ω<sup>p</sup>* is the plasma frequency of the conductor (depending on the dielectric constant in vacuum, the density of free electrons, electron mass, and charge), *ω* is the angular frequency of the incidence wave, and *γ* is the damping coefficient of the conductor. The plasma frequency for aluminum is *<sup>ω</sup><sup>p</sup>* = 2.123 <sup>×</sup> <sup>10</sup><sup>16</sup> rad/s, while the damping coefficient is *<sup>γ</sup>* = 1.975 <sup>×</sup> <sup>10</sup>11/s. To reduce the computational complexity of the developed model (instead of the simulation of a large matrix of structural elements) just a single-element analysis was performed. This simplification was based on the application of Floquet periodic boundary conditions on the side boundaries. Such simplification of the model is acceptable, as will be shown later in this paper, through the comparison of modeled and measured resonant frequencies and *S*<sup>21</sup> reflection coefficient frequency responses.

respectively [25]. The Drude model was used to describe the value of the aluminum re-

on the bottom boundary of the computational model in order to fully absorb the outgoing wave and omit its influence on the results including the Fabry–Perot effect. However, this **Figure 4.** Numerical model of SRR metasurface on the silicon wafer (the thickness of the Si substrate not in scale). **Figure 4.** Numerical model of SRR metasurface on the silicon wafer (the thickness of the Si substrate not in scale).

effect can also occur in the case of the Si substrate. For this reason, the solution proposed in paper [27] was utilized: the Si wafer consisting of metasurface (not in scale in Figure 4) was backed with an Si layer (not an air layer such as in ordinary operating conditions), and the receiving port was placed on the back surface of the Si wafer. This causes a lack of reflection on the Si wafer's back layer boundary. Thus, this configuration makes it possible to eliminate the Fabry–Perot effect while conserving the original response of the metasurface. The Fabry–Perot effect strongly affects the results obtained in the terahertz To investigate the influence of the fabrication process and material properties on the quality of the designed SRR structure, the following parameters were analyzed: etching deviation *e*, metallization thickness *t*Al, and SiO2 layer thickness *t*SiO2. The quality of the metasurface affected by the above parameters was determined by the deviation of resonant frequency Δ*f*r: ∆୰ ൌ ୰ୈ െ ୰ (3) The wave excitation and scattering matrix parameters calculation were enabled through ports (excitation and *S*<sup>11</sup> calculation on port 1, *S*<sup>21</sup> calculation on port 1 and 2). Port 1 is placed at some distance over the MS element separated by an air layer, while port 2 is placed on the end of the metasurface stack. The scattering boundary condition was set on the bottom boundary of the computational model in order to fully absorb the outgoing wave and omit its influence on the results including the Fabry–Perot effect. How-

where *f*rD—resonant frequency of the structure with deviated parameters, *f*rB—resonant

The developed numerical model was utilized to perform the considered analyses. The exemplary results of the simulations are presented in Figure 5; the frequency response of transmission |*S*21| and the reflection |*S*11| scattering matrix parameters were presented for the dimensions given in Table 1. One can observe that there are two resonances: 1 THz and 2.77 THz. The distribution of the electric field on the surface of MS is shown for the

determined using a developed finite element method (FEM)-based numerical model.

**3. Results**

ever, this effect can also occur in the case of the Si substrate. For this reason, the solution proposed in paper [27] was utilized: the Si wafer consisting of metasurface (not in scale in Figure 4) was backed with an Si layer (not an air layer such as in ordinary operating conditions), and the receiving port was placed on the back surface of the Si wafer. This causes a lack of reflection on the Si wafer's back layer boundary. Thus, this configuration makes it possible to eliminate the Fabry–Perot effect while conserving the original response of the metasurface. The Fabry–Perot effect strongly affects the results obtained in the terahertz frequency range (multiple, relatively high amplitude ripples obtained in spectrum); thus, in the case of frequency domain computations, some sophisticated methods for its reduction should be applied.

In this high-frequency electromagnetic problem the following frequency domain equation is solved [28,29]:

$$\nabla \times \mu\_{\rm r}^{-1}(\nabla \times E) - k\_0^2 \left( \varepsilon\_{\rm r} - \frac{j\sigma\_{\varepsilon}}{\omega \varepsilon\_0} \right) E = 0 \tag{2}$$

where: *E*—electric field, *µ*r—relative permeability, *ε*0—permittivity of vacuum, *ε*r—relative permittivity, *σe*—electric conductivity, *k*0—wave number, *ω*—angular frequency. The exemplary mesh used during numerical calculations consists of 106,029 domain elements, 17,278 boundary elements, and 1602 edge elements.

To reduce the computational complexity of the developed model (instead of the simulation of a large matrix of structural elements) just a single-element analysis was performed. This simplification was based on the application of Floquet periodic boundary conditions on the side boundaries. Such simplification of the model is acceptable, as will be shown later in this paper, through the comparison of modeled and measured resonant frequencies and *S*<sup>21</sup> reflection coefficient frequency responses.

To investigate the influence of the fabrication process and material properties on the quality of the designed SRR structure, the following parameters were analyzed: etching deviation *e*, metallization thickness *t*Al, and SiO<sup>2</sup> layer thickness *t*SiO2. The quality of the metasurface affected by the above parameters was determined by the deviation of resonant frequency ∆*f*r:

$$
\Delta f\_{\mathbf{r}} = f\_{\mathbf{r}\mathbf{D}} - f\_{\mathbf{r}\mathbf{B}} \tag{3}
$$

where *f*rD—resonant frequency of the structure with deviated parameters, *f*rB—resonant frequency of the structure with a base set of parameters (shown in Table 1).

The influence of each of the above parameters on the behavior of the metasurface was determined using a developed finite element method (FEM)-based numerical model.

#### **3. Results**

The developed numerical model was utilized to perform the considered analyses. The exemplary results of the simulations are presented in Figure 5; the frequency response of transmission |*S*21| and the reflection |*S*11| scattering matrix parameters were presented for the dimensions given in Table 1. One can observe that there are two resonances: 1 THz and 2.77 THz. The distribution of the electric field on the surface of MS is shown for the first resonance (frequency of 1 THz). As was expected, there is an enhancement of the electric field intensity in the vicinity of the capacitive gap. The presented results show that the designed metasurface meets the requirements assumed in the design process. This positive verification allowed us to perform the fabrication process based on the dimensions presented in Table 1. The photo of the fabricated structure is shown in Figure 6.

first resonance (frequency of 1 THz). As was expected, there is an enhancement of the electric field intensity in the vicinity of the capacitive gap. The presented results show that the designed metasurface meets the requirements assumed in the design process. This positive verification allowed us to perform the fabrication process based on the dimensions presented in Table 1. The photo of the fabricated structure is shown in Figure 6.

first resonance (frequency of 1 THz). As was expected, there is an enhancement of the electric field intensity in the vicinity of the capacitive gap. The presented results show that the designed metasurface meets the requirements assumed in the design process. This positive verification allowed us to perform the fabrication process based on the dimensions presented in Table 1. The photo of the fabricated structure is shown in Figure 6.

*Electronics* **2023**, *12*, x FOR PEER REVIEW 6 of 13

**Figure 5.** Results of FEM numerical simulations: frequency response (**left**) and electric field distribution on the surface of SRR (**right**). **Figure 5.** Results of FEM numerical simulations: frequency response (**left**) and electric field distribution on the surface of SRR (**right**). **Figure 5.** Results of FEM numerical simulations: frequency response (**left**) and electric field distribution on the surface of SRR (**right**).

**Figure 6.** Photo of manufactured SRR metasurface on the silicon wafer (light area—Al metallization; dark area—lack of metallization, SiO2). **Figure 6.** Photo of manufactured SRR metasurface on the silicon wafer (light area—Al metallization; dark area—lack of metallization, SiO2). **Figure 6.** Photo of manufactured SRR metasurface on the silicon wafer (light area—Al metallization; dark area—lack of metallization, SiO<sup>2</sup> ).

The fabricated structure was experimentally evaluated to verify the results of the design process, as well as the numerical model and the applied assumptions (such as frequency domain computation, Floquet periodic boundary conditions, plane wave excitation, dispersive modeling of the metallization layer, etc.). The measurements were performed using the THz-TDS (terahertz time domain spectroscopy) system TeraFlash smart of Toptica, Germany, shown in Figure 7. A transmission setup was utilized, whereby the transmission Tx and receiver Rx heads are on opposite sides of the examined device/structure under test (DUT). The electromagnetic pulse with a terahertz frequency content is generated by a photoconductive antenna and focused on the evaluated metasurface structure using a set of curved reflectors. After interaction with the DUT, the modified pulse defocuses and propagates through another set of curved reflectors and is passed to the receiving photoconductive antenna. The focusing spot size (the area of incidence of the The fabricated structure was experimentally evaluated to verify the results of the design process, as well as the numerical model and the applied assumptions (such as frequency domain computation, Floquet periodic boundary conditions, plane wave excitation, dispersive modeling of the metallization layer, etc.). The measurements were performed using the THz-TDS (terahertz time domain spectroscopy) system TeraFlash smart of Toptica, Germany, shown in Figure 7. A transmission setup was utilized, whereby the transmission Tx and receiver Rx heads are on opposite sides of the examined device/structure under test (DUT). The electromagnetic pulse with a terahertz frequency content is generated by a photoconductive antenna and focused on the evaluated metasurface structure using a set of curved reflectors. After interaction with the DUT, the modified pulse defocuses and propagates through another set of curved reflectors and is passed to the receiving photoconductive antenna. The focusing spot size (the area of incidence of the The fabricated structure was experimentally evaluated to verify the results of the design process, as well as the numerical model and the applied assumptions (such as frequency domain computation, Floquet periodic boundary conditions, plane wave excitation, dispersive modeling of the metallization layer, etc.). The measurements were performed using the THz-TDS (terahertz time domain spectroscopy) system TeraFlash smart of Toptica, Germany, shown in Figure 7. A transmission setup was utilized, whereby the transmission Tx and receiver Rx heads are on opposite sides of the examined device/structure under test (DUT). The electromagnetic pulse with a terahertz frequency content is generated by a photoconductive antenna and focused on the evaluated metasurface structure using a set of curved reflectors. After interaction with the DUT, the modified pulse defocuses and propagates through another set of curved reflectors and is passed to the receiving photoconductive antenna. The focusing spot size (the area of incidence of the terahertz beam) is less than 1 mm diameter. The results of the measurements are presented in Figure 8. The presented time domain pulses were obtained for two scenarios: the Si wafer without metasurface (reference signal), and the Si wafer containing an MS. As the frequency domain data are of more importance in the considered analysis (resonant frequency analysis), a transformation of the data from the time to frequency domain was obtained using Fourier transform. In the first step, the measurement without DUT was made (air transmission

spectrum) in order to verify the transmission band. In the obtained frequency response, the usable signal (above the noise level) can be observed for the frequency band below 2.5 THz. This means that the second resonance will not be possible to measure for this configuration (the frequency responses presented in Figure 8 obtained for the Si wafer also confirm a maximum usable frequency of 2.5 THz). Within the transmission band, there are multiple absorption peaks visible, caused by the absorption of certain frequencies by water, nitrogen, and oxygen molecules. In the next step, the measurement of the Si wafer without metallization was made to obtain the reference signal. Finally, the measurement of full DUT (the Si wafer with a metallization metasurface) was performed. One can observe the decrease in the transmission spectrum in the vicinity of a frequency of 1 THz, which is caused by the resonant behavior of MS. After de-embedding, the metasurface transmission spectrum consists of just one (primary) resonance for the desired frequency of 1 THz. The resonant frequency obtained in the simulations was *f*r1Sim = 1.0023 THz, and in the measurements was *f*r1Meas = 1.040 THz. Good agreement between the simulations and measurement results was proved by an absolute difference of 37.7 GHz and a relative difference of 3.76%. The measurement results are consistent with the resonance frequency *f*r = 1 THz assumed in the design process and is within the deviation of +46.47 GHz/−42.83 GHz determined (in a later part of this work) by numerical analysis. As indicated in the article, due to the limited bandwidth of the measurement system used, experimental verification in the case of the second resonance was not possible. made (air transmission spectrum) in order to verify the transmission band. In the obtained frequency response, the usable signal (above the noise level) can be observed for the frequency band below 2.5 THz. This means that the second resonance will not be possible to measure for this configuration (the frequency responses presented in Figure 8 obtained for the Si wafer also confirm a maximum usable frequency of 2.5 THz). Within the transmission band, there are multiple absorption peaks visible, caused by the absorption of certain frequencies by water, nitrogen, and oxygen molecules. In the next step, the measurement of the Si wafer without metallization was made to obtain the reference signal. Finally, the measurement of full DUT (the Si wafer with a metallization metasurface) was performed. One can observe the decrease in the transmission spectrum in the vicinity of a frequency of 1 THz, which is caused by the resonant behavior of MS. After de-embedding, the metasurface transmission spectrum consists of just one (primary) resonance for the desired frequency of 1 THz. The resonant frequency obtained in the simulations was *f*r1Sim = 1.0023 THz, and in the measurements was *f*r1Meas = 1.040 THz. Good agreement between the simulations and measurement results was proved by an absolute difference of 37.7 GHz and a relative difference of 3.76%. The measurement results are consistent with the resonance frequency *f*r = 1 THz assumed in the design process and is within the deviation of +46.47 GHz/−42.83 GHz determined (in a later part of this work) by numerical analysis. As indicated in the article, due to the limited bandwidth of the measurement system used, experimental verification in the case of the second resonance was not possible.

terahertz beam) is less than 1 mm diameter. The results of the measurements are presented in Figure 8. The presented time domain pulses were obtained for two scenarios: the Si wafer without metasurface (reference signal), and the Si wafer containing an MS. As the frequency domain data are of more importance in the considered analysis (resonant frequency analysis), a transformation of the data from the time to frequency domain was obtained using Fourier transform. In the first step, the measurement without DUT was

*Electronics* **2023**, *12*, x FOR PEER REVIEW 7 of 13

**Figure Figure 7. 7.** Photo Photo of THz TDS measuring system. of THz TDS measuring system.

After positive verification of the developed numerical model and all the assumptions made, analysis of the influence of selected geometrical parameters was performed. These parameters (etching deviation *e*, metallization thickness *t*Al, and SiO<sup>2</sup> layer thickness *t*SiO2) were selected as the most susceptible to possible deviations during the production process. The values of the *t*Al and *t*SiO2 parameters presented in Table 1 were taken as reference values, and during the analysis, these values were changed in the range of ±20% of the reference value: *e* = −20 . . . + 20%, *t*Al = 0.400 . . . 0.600 µm and *t*SiO2 = 0.136 . . . 0.204 µm. The etching deviation parameter is a more complicated value, consisting of the line width *w*, which was changed from the range of ±20%, while simultaneously changing the gap *g* and the square loop size *a*. The considered values of *w*, *g*, and *a* for the different values of the etching deviation parameter *e* are shown in Table 2.

**Figure 8.** Results of measurements using THz TDS spectroscope. Results were obtained through the Fourier transform of time domain data. **Figure 8.** Results of measurements using THz TDS spectroscope. Results were obtained through the Fourier transform of time domain data.


After positive verification of the developed numerical model and all the assumptions **Table 2.** List of parameter sets for individual etching deviation parameter values *e*.

**Etching Deviation**  *e* **[%] Line Width**  *w* **[µm] Gap**  *g* **[µm] Square Loop Size**  *a* **[µm]** −20 2 3 21.2 −15 2.125 2.875 21.325 −10 2.25 2.75 21.45 An exemplary result of such analysis—the shift of the resonant curve caused only by the change in the metallization layer thickness *t*Al—is presented in Figure 9. One can observe the increase in resonant frequency *f*r1 while the metallization thickness *t*Al increases. This relationship is presented for both resonances and all three considered parameters in Figure 10. In all cases we can observe a quasi-linear relationship when considering ±20% variation ranges.

−5 2.375 2.625 21.575 0 2.5 2.5 21.7 5 2.625 2.375 21.825

**Table 2.** List of parameter sets for individual etching deviation parameter values *e*.

**Figure 9.** Results of FEM numerical simulations: shift of resonant curve caused by the thickness of the metallization layer *t*Al. **Figure 9.** Results of FEM numerical simulations: shift of resonant curve caused by the thickness of the metallization layer *t*Al. **Figure 9.** Results of FEM numerical simulations: shift of resonant curve caused by the thickness of the metallization layer *t*Al.

**Figure 10.** Numerically calculated influence of various geometrical parameters (*e*, *t*Al, and *t*SiO2) on the shift of resonant frequencies Δ*f*r1 and Δ*f*r2: etching *e* (**left**), aluminum metallization thickness *t*Al (**right**), and silicon dioxide layer thickness *t*SiO2 (**bottom**). **Figure 10.** Numerically calculated influence of various geometrical parameters (*e*, *t*Al, and *t*SiO2) on the shift of resonant frequencies Δ*f*r1 and Δ*f*r2: etching *e* (**left**), aluminum metallization thickness *t*Al (**right**), and silicon dioxide layer thickness *t*SiO2 (**bottom**). **Figure 10.** Numerically calculated influence of various geometrical parameters (*e*, *t*Al, and *t*SiO2) on the shift of resonant frequencies ∆*f*r1 and ∆*f*r2: etching *e* (**left**), aluminum metallization thickness *t*Al (**right**), and silicon dioxide layer thickness *t*SiO2 (**bottom**).

**Table 3.** Comparison of various parameters' influences on the maximum shift of resulting resonant frequencies. **Uncertainty Source Δ***f***r1 [GHz] Δ***f***r2 [GHz] Table 3.** Comparison of various parameters' influences on the maximum shift of resulting resonant frequencies. **Uncertainty Source Δ***f***r1 [GHz] Δ***f***r2 [GHz]** The influence of the analyzed parameters can be considered based on an equivalent circuit—for SRR—schematically modeled as a parallel resonant circuit with a resonance frequency *f*rSCH [30]:

$$f\_{\rm rSCH} = \frac{1}{2\pi\sqrt{LC}}\tag{4}$$

Aluminum thickness, *t*Al® − 20% −10.02 −23.48 Aluminum thickness, *t*Al® + 20% 9.82 24.80 Aluminum thickness, *t*Al® − 20% −10.02 −23.48 Aluminum thickness, *t*Al® + 20% 9.82 24.80 where *f*rSCH—resonant frequency of equivalent circuit, *L*—equivalent inductance,*C*—equivalent capacitance.

Silicon dioxide thickness, *t*SiO2® − 20% −12.76 −35.91 Silicon dioxide thickness, *t*SiO2® + 20% 12.98 36.02 First worst case (*e*) – 20%, (*t*Al + *t*SiO2) + 20% 46.47 170.61 Silicon dioxide thickness, *t*SiO2® − 20% −12.76 −35.91 Silicon dioxide thickness, *t*SiO2® + 20% 12.98 36.02 Increasing the thickness of the metallization increases the conductance but decreases the inductance of the wire constituting the resonant ring. Reducing the inductance in the equivalent circuit increases the resonant frequency.

Second worst case (*e*) + 20%, (*t*Al + *t*SiO2) − 20% −42.83 −192.12 First worst case (*e*) – 20%, (*t*Al + *t*SiO2) + 20% 46.47 170.61 Second worst case (*e*) + 20%, (*t*Al + *t*SiO2) − 20% −42.83 −192.12 The capacitance of the SRR structure depends to the greatest extent on the capacitance of the gap, which can be approximated by a flat capacitor, in which the electric field is partially distributed in air and partially in a silicon wafer, in particular the SiO<sup>2</sup> layer on the surface. Therefore, to determine the capacitance, it is necessary to use the equivalent permittivity (depending on the depth of the SiO<sup>2</sup> layer). In the case of a thin SiO<sup>2</sup> layer, part of the electric field exists in the Si (high permeability) layer. The increase in the thickness of the SiO<sup>2</sup> layer causes the electric field to be distributed to a greater extent in the material characterized by lower permittivity, so as a consequence the resultant capacitance decreases. A smaller value of the gap capacitance increases the resonant frequency.

According to Table 2, an increase in the etching deviation factor causes an increase in linewidth and square loop size and a decrease in capacitive gap size. An increasing linewidth causes a decrease in the inductance and, as a result, an increase in resonant frequency. A reduction in gap size causes a higher capacitance and consequently a lower resonant frequency. If the square loop size increases, the inductance increases too, thus the resonant frequency decreases. Finally, in the case of the etching deviation factor increasing, the resulting resonance moves to lower values.

The comparison of the influence of considered parameters on the maximum shift of resulting resonant frequencies is presented in Table 3. The comparison was divided into two stages: first, each of the considered parameters was analyzed separately, then two worst-case scenarios were verified:


**Table 3.** Comparison of various parameters' influences on the maximum shift of resulting resonant frequencies.


This was possible because of similar and monotonic characteristics of the relationships shown in Figure 10.

In the case of only one parameter variation, the biggest influence was caused by the line etching deviation *e*, while the smallest one was by aluminum thickness *t*Al. The frequency shifts, as expected, were 2–5 times higher for the second resonance. Finally, in the case of the worst-case scenario (all the considered parameters deviations were taken into account), the maximum change in *f*r1 was close to 46 GHz and the change in *f*r2 was 192 GHz. This caused the following relative deviations of resonant frequencies: 4.6% for *f*r1 and 6.9% for *f*r2. The shifts of resonant curves caused by the cumulative change in etching deviation *e*, aluminum layer thickness *t*Al, and silicon dioxide thickness *t*SiO2 are presented in Figure 11.

In the case of only one parameter variation, the biggest influence was caused by the line etching deviation *e*, while the smallest one was by aluminum thickness *t*Al. The frequency shifts, as expected, were 2–5 times higher for the second resonance. Finally, in the case of the worst-case scenario (all the considered parameters deviations were taken into account), the maximum change in *f*r1 was close to 46 GHz and the change in *f*r2 was 192 GHz. This caused the following relative deviations of resonant frequencies: 4.6% for *f*r1 and 6.9% for *f*r2. The shifts of resonant curves caused by the cumulative change in etching deviation *e*, aluminum layer thickness *t*Al, and silicon dioxide thickness *t*SiO2 are presented

**Figure 11.** Results of FEM numerical simulations: shift of resonant curve caused by the cumulative change in etching deviation *e*, aluminum layer thickness *t*Al, and silicon dioxide thickness *t*SiO2 worst cases. **Figure 11.** Results of FEM numerical simulations: shift of resonant curve caused by the cumulative change in etching deviation *e*, aluminum layer thickness *t*Al, and silicon dioxide thickness *t*SiO2—worst cases.

#### **4. Conclusions 4. Conclusions**

in Figure 11.

In this paper, the design and fabrication of split-ring-resonator-based metasurfaces were performed and a further analysis of the fabrication process uncertainty on resulting MS parameters was numerically conducted. The process of terahertz metasurface fabrication is relatively complicated (lots of parameters need to be precisely controlled), and may be the source of some deviations from the planned (designed) dimensions resulting in shifts of assumed electromagnetic parameters. The sources of probable inaccuracies in the fabrication process were identified (etching process deviation *e*, aluminum metallization thickness *t*Al, and silicon dioxide layer thickness *t*SiO2 variations) and their influence was analyzed using the numerical model. The commonly utilized metasurface structure of the square split-ring resonator was chosen and designed for operation in the terahertz frequency range. In this case, not the full matrix of structural elements but a single-element numerical analysis (with periodic boundary conditions) was utilized. This assumption regarding the periodicity enables the examination of various fabrication aspects, as long as there is no asymmetry in the metasurface's geometry. In this paper, the design and fabrication of split-ring-resonator-based metasurfaces were performed and a further analysis of the fabrication process uncertainty on resulting MS parameters was numerically conducted. The process of terahertz metasurface fabrication is relatively complicated (lots of parameters need to be precisely controlled), and may be the source of some deviations from the planned (designed) dimensions resulting in shifts of assumed electromagnetic parameters. The sources of probable inaccuracies in the fabrication process were identified (etching process deviation *e*, aluminum metallization thickness *t*Al, and silicon dioxide layer thickness *t*SiO2 variations) and their influence was analyzed using the numerical model. The commonly utilized metasurface structure of the square split-ring resonator was chosen and designed for operation in the terahertz frequency range. In this case, not the full matrix of structural elements but a single-element numerical analysis (with periodic boundary conditions) was utilized. This assumption regarding the periodicity enables the examination of various fabrication aspects, as long as there is no asymmetry in the metasurface's geometry.

It was assumed that the geometrical parameters (*e*, *t*Al, and *t*SiO2) were modified in the range of ± 20%. All the evaluated sources of geometric deviation caused a linear change in both resonant frequencies. This caused the following relative deviations of resonant frequencies: 4.6% for *f*r1 and 6.9% for *f*r2. Such frequency deviations can cause significant degradation of the MS performance in some applications. Depending on the application of a particular metasurface, other parameters (than resonant frequency) should be considered for evaluation; e.g., in the case of reconfigurable structures (flexible or MEMS devices), It was assumed that the geometrical parameters (*e*, *t*Al, and *t*SiO2) were modified in the range of ±20%. All the evaluated sources of geometric deviation caused a linear change in both resonant frequencies. This caused the following relative deviations of resonant frequencies: 4.6% for *f*r1 and 6.9% for *f*r2. Such frequency deviations can cause significant degradation of the MS performance in some applications. Depending on the application of a particular metasurface, other parameters (than resonant frequency) should be considered for evaluation; e.g., in the case of reconfigurable structures (flexible or MEMS devices), tunability range can be evaluated.

**Author Contributions:** Conceptualization, P.L., M.H., U.M. and A.K.; methodology, P.L., M.H., U.M. and A.K.; process development and fabrication, U.M., A.K. and A.F.; numerical modelling, M.H., P.G. and P.L.; measurements, P.L. and M.H.; writing—original draft preparation, P.L, M.H. and A.K.; writing—review and editing, P.L. and A.K.; project administration, P.L. and U.M.; fund-**Author Contributions:** Conceptualization, P.L., M.H., U.M. and A.K.; methodology, P.L., M.H., U.M. and A.K.; process development and fabrication, U.M., A.K. and A.F.; numerical modelling, M.H., P.G. and P.L.; measurements, P.L. and M.H.; writing—original draft preparation, P.L, M.H. and A.K.; writing—review and editing, P.L. and A.K.; project administration, P.L. and U.M.; funding acquisition, P.L. and U.M. All authors have read and agreed to the published version of the manuscript.

ing acquisition, P.L. and U.M. All authors have read and agreed to the published version of the manuscript. **Funding:** This work was conducted within the project "RE-TERA- Reconfigurable terahertz devices for EM waves manipulation and sensing applications", co-financed by the Polish National Agency for Academic Exchange (NAWA, Poland) and German Academic Exchange Service (DAAD, Germany) under the grant no. PPN/BDE/2021/1/00012/U/00001/ID: 57602825.

**Acknowledgments:** We would like to thank Xenia Seng from the Institute of Microsystems Technology, Furtwangen University, for carrying out the microfabrication process.

**Conflicts of Interest:** The authors declare no conflict of interest.

tunability range can be evaluated.

## **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
